Investigations of quasi-static vortex structures in 2D sand specimen under passive earth pressure conditions based on DEM and Helmholtz–Hodge vector field decomposition

The paper presents some two-dimensional simulation results of granular vortex-structures in cohesionless initially dense sand during a quasi-static passive wall translation. 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. In order to detect vortex-structures, the Helmholtz–Hodge decomposition of a vector field from DEM calculations was used. This approach enabled us to distinguish both incompressibility and vorticity in the granular displacement field. In addition the predominant periods of vortices during horizontal wall movement were determined. The vortices were strongly connected to shear localization. They localized in locations where shear zones ultimately developed. In addition, the vortex-structures were calculated during plane strain compression.


Introduction
Granular vortex-structures (swirling motion of several grains around its central point) were frequently observed in experiments on granular materials [1][2][3][4] and in calculations using the discrete element method (DEM) [5][6][7][8][9][10][11][12][13][14][15][16]. The experi-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 ments were performed during Couette shearing [1], plane strain compression [2] and simple shear [3,4]. The vortices were observed to form in association with the onset of peak stress. They appeared only occasionally, quickly dissipating [1][2][3][4]. The granular vortex-structures were also obtained with the aid of the finite element method [17]. They became apparent in experiments [1][2][3][4] and calculations (e.g. [5,[13][14][15][16]) when the motion associated with uniform (affine) strain was subtracted from the actual granular deformation. They are reminiscent of turbulence in fluid dynamics [5], however the amount of the grain rotation is several ranges of magnitude smaller (∼0.01 • -0.1 • ) than the fluid vortex rotation and granular flow is too slow to induce inertial forces characteristic for turbulences in fluid. According to Peters and Walizer [13] vortices represent an independent flow field following its own governing equations and satisfying its own (null) boundary conditions. Tordesillas et al. [11] showed that two classes of vortices emerging: primary ones concentrated in the shear zone and secondary ones forming next to the zone boundaries. A dominant mechanism responsible for the vortex formation was the breakage of force chains [11,16]. The collapse of main force chains lead to a formation of larger voids and their build-up to a formation of smaller voids [11,16]. Vortex dynamics were consistent with stick-slip dynamics [16]. The vortices have been mainly observed in shear zones [11,16,18] which are the fundamental phenomenon in granular bodies. A precise mechanism behind the vortex evolution with its connection to shear localization in granular bodies still remains elusive. In numerical calculations, shear localization is usually identified in granular bodies by grain rotations (in DEM) or micro-polar rotations or by an increase of void ratio (in FEM) [19].
The objective of this paper is to report the results of comprehensive 2D studies by DEM on vortex-structures in sand behind a rigid wall during its quasi-static passive translation by using the Helmholtz-Hodge decomposition (HHD) of a vector field [20,21]. Attention was paid to the relationship between vortex-structures and shear localization with respect to the location and formation moment. The analyses were carried out with spheres with contact moments [16] to approximately capture the irregular grain shape. 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. A three-dimensional discrete model YADE developed at University of Grenoble by F.V. Donze and his co-workers was applied [22][23][24][25][26][27]. The discrete calculations were solely carried out with initially dense sand.
In our previous paper we calculated 2D vortex-and anti-vortex structures during a quasi-static plane strain compression test solved by the DEM [28]. Two-dimensional vortex-and anti-vortex structures were determined by a method based on orientation angles of displacement fluctuation vectors of neighbouring single spheres [28]. The proposed method used for the detection of 2D granular vortex/anti-vortex-structures was very effective. The method detected all vortex and anti-vortex-structures regardless of the displacement vector length. The 2D vortex-and antivortex structures solely appeared in the main inclined shear zone. Their occurrence did not depend on the specimen depth. The anti-vortices turned out be the best precursor of the location of shear zones since they appeared from the beginning of the deformation process, i.e. significantly earlier than e.g. based on the average cumulative grain rotation. Other results showed that the increase of the grain-non-regularity decreased the predominant period of left-handed vortices and anti-vortices by 25-40%. This method had however 3 disadvantages: (1) it operated on displacement fluctuation vectors (not directly on displacement increment vectors), (2) it was solely designed for 2D granular flow and (3) it depended on several parameters such as: region size of the average background translation, averaging region size of local displacement increments, searching circle radius and grid move distance [18]. The method used in this paper has the following advantages over the previous one: (1) it operates directly on displacement increment vectors, (2) it is designed for identifying vortices in both two-and three-dimensional kinematic fields, (3) it does not need any additional parameters for calculations and (4) it can be used for DEM and FEM calculations. The method does not determine the size of vortex-structures. There exist also the methods which take the vortex size into account [29]. However this size depends on the displacement vector lengths, displacement vector rotations, ratio between the height and width of vectors assumed for describing the vortex shape and number of grains comprising the vortex. In the paper [16] in order to mathematically describe 2D vortices, the displacement fluctuation vector of each grain in the neighbourhood of each central grain was decomposed into 2 vectors: the normal and tangential to the movement direction of each central grain. Only the tangential displacement fluctuations were assumed to be responsible for a vortex, i.e. if the neighbouring grains had solely the tangential displacement fluctuation component, the central grain was assumed to be located in a vortex mid-point. In order to avoid the prescribing pure shearing along a shear zone as a vortex, a certain limitation was imposed, i.e. the condition that the sum of the tangential displacement fluctuations had to be at least twice as big as the sum of the normal displacement fluctuations (to be classified as a part of the vortex). It was assumed in Peters and Walizer [13] that a 2D vortex was connected with the rigid body rotation (spin) in the displacement fluctuation field. During 2D calculations described in Tordesillas et al. [11], vortex cores were identified by partitioning a vector displacement field into triangles, each analyzed with respect to a direction-spanning property (the vectors at the triangle vertices were located within a certain range of angles, i.e. each vector pointed to a unique direction range). A vortex was then identified for each vortex core with the vortex boundary set by the last member vector rotating in the vortex direction. The full extent of each vortex was established by a systematic check of the direction of neighbouring vectors (in the second ring, third ring, etc.) at increasing radial distance from the vortex core until a vector was found that rotated in an opposite direction to the corresponding vortex.
Our paper consists of two main parts. In the first part, discrete elements results of a passive earth pressure problem were shortly described to show the capability of DEM to realistically simulate the pattern of shear zones in sand [17]. In the second part, the formation of 2D vortex-structures by the HHD was discussed based on direct grain displacements. In addition, the vortex-structures were calculated during plane strain compression.

Three-dimensional DEM model
In order to simulate the behaviour of real sand, the 3D spherical discrete model YADE, developed at University of Grenoble [22][23][24][25][26][27], was used by taking 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. A choice of a very simple linear elastic normal contact was intended to capture on average various contacts possible in real sands. The normal and tangential forces were linked to the displacements through the normal stiffness K n and tangential stiffness K s (Fig. 1A) Fig. 1 Mechanical response of linear contact model without (A) and with contact moments (A+B) [22,23]: (a) tangential contact model, (b) 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, X s -tangential displacement vector, ω-angular rotation vector, μ-inter-particle friction angle, η-limit rolling coefficient where U is the penetration depth between discrete elements, N denotes the unit normal vector at the contact point and X s is the incremental tangential displacement vector. The unloading was assumed to be purely elastic (Fig. 1C). The stiffness parameters were calculated in terms of the modulus of elasticity of the grain contact E c and two contacting grain radii R A and R B (to determine the normal stiffness K n ) and in terms of the modulus of elasticity E c and Poisson's ratio υ c of the grain contact, and grain radii R A and R B (to determine the tangential stiffness K s ) of two contacting spheres, respectively [22] If the grain radius R A = R B = R, the stiffness parameters are equal to: K n = E c R and K s = υ c E c R (thus K s /K n = υ c ), respectively. The frictional sliding starts at the contact point when the contact forces F s and F n satisfy the limit Coulomb condition (Fig. 1a) with μ as the inter-particle friction angle (tension was not allowed). No forces are transmitted when grains are 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 ν [30,31]. In order to increase the rolling resistance of pure spheres, clusters of spheres or contact moments were introduced [30]. The normal force was assumed to contribute to the rolling resistance. The contact moment increments were calculated by means of the rolling stiffness K r multiplied by the angular rotational increment vectors ω (Fig. 1B) In Eq. 5, the angular rotational increment vectors do not depend on the spherical grain radii in contrast to equations which take the different radii into account [32,33].
The rolling stiffness K r [kNm] in Eq. 5 was related to the tangential stiffness K s [kN/m] in Eq. 2 by the following formula proposed by Iwashita and Oda [34] where β is the dimensionless rolling stiffness coefficient and R is the equivalent grain radius (at small displacements dX r ≈ dX s ). The dimensionless rolling coefficient η specifies the limit friction moment of the rolling motion In order to dissipate excessive kinetic energy in the discrete system, a simple local non-viscous damping scheme was adopted [35], by assuming a change of forces and moment reduced due to the damping effect specified by the parameter α and where F k and M k are the kth components of the residual force and moment vector and ν k and ω k are the kth components of the translational and rotational velocity. A positive damping coefficient α is smaller than 1 (sgn(•) returns the sign of the kth component of velocity). The equations are separately applied to each k-th component of a 3D vector x, y and z.
Note that the effect of damping is insignificant in quasi-static calculations [31]. Although a non-linear contact law is more realistic, a linear contact law provides similar results with the significantly reduced computation time [31] and therefore was used in the present simulations. The five main local material parameters are necessary in our DEM simulations: E c (modulus of elasticity of the grain contact), ν 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 , ν 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 [36]. The procedure for determining the material parameters in DEM was described in detail by Kozicki et al. [30,31]. The index properties of Karlsruhe sand are: mean grain diameter d 50 = 0.50 mm, grain size between 0.08 and 1.8 mm, uniformity coefficient U c = 2, maximum specific weight γ max d = 17.4 kN/m 3 , minimum void ratio e min = 0.53, minimum specific weight γ min d = 14.6 kN/m 3 and maximum void ratio e max = 0.84. 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. The constants E c and ν c do not correspond to the elastic constants of grains [31]. Note that the other set of the material constants μ, β and η is also possible.

DEM results of passive earth pressure model tests
The DEM calculation results were described in detail in Nitka et al. [16]. The simulations were performed for a 2D sand body of l w = 0.40 m length and h w = 0.20 m height in order to compare with experiments with Karlsruhe sand (d 50 = 0.5 mm) [37,38]. The vertical retaining wall and the bottom of the granular specimen were assumed to be stiff and very rough, i.e. there was no relative displacement along vertical and bottom surface [19]. Since the experiments were idealized as a 2D boundary value problem and the effect of the specimen depth in the out of plane direction turned out to be almost negligible during direct shearing in DEM calculations [15] in order to significantly accelerate simulations, the computations were performed with the specimen depth equal to the grain size (i.e. one layer of spheres was simulated along the depth only).
The spheres with d 50 = 1.0 mm, characterized by a linear grain size distribution, were assumed (grain size range 0.5-1.5 mm, 62,600 spheres). Along the height of 200 mm, about 200 spheres were used. The initial void ratio of sand, obtained by generating random spheres above a box and then allowing them to fall down by gravity, was e o = 0.62. The loading speed was slow enough to ensure that the tests were conducted under quasi-static conditions. The calculated mean inertial number (which quantifies the significance of dynamic effects) in the mid-length of the curved shear zone for the maximum horizontal earth pressure force was in the analyses I =γ d 50 P ρ = 5 · 10 −4 (γ = 0.751/s-the shear rate, P = 5 kPa-the mean pressure and ρ = 2.55 g/cm 3 ). The inertial number obviously changed along the specimen height and in time during the entire granular flow. The value of I ≤ 10 −3 corresponds usually to a quasi-static regime [39]. Figure 2A shows the evolution of the resultant normalized horizontal earth pressure force (earth pressure coefficient) wall displacement u h / h w (h w = 0.2 m, E h -the horizontal force acting on the wall, γ = 16.75 kN/m 2 -the initial unit weight) for d 50 = 1 mm from DEM simulations. The normalized horizontal earth pressure force evolved typically for initially dense granulates in biaxial compression, triaxial compression and direct shearing. The specimen exhibited the initial strain hardening up to the peak (u h / h w = 0.038), followed by some softening before the common asymptote was reached (critical state). The horizontal force fluctuated after the peak that was attributed to the build-up and collapse of force chains-the main carrier of stresses transferred within the granular assembly [16]. The earth pressure coefficient was K max p = 30 for d 50 = 1 mm. It can be thus anticipated that for d 50 = 0.5 mm (real sand), K max p should be about 25-27. The value of K max p = 30 for d 50 = 1 mm was slightly smaller than K max p = 31 obtained by FEM (d 50 = 0.5 mm) [38] and was closer to the engineering earth pressure coefficients [40].
The distribution of single sphere rotations ω during wall translation is presented in Fig. 2B which are usually the best indicator for shear localization in DEM (red denotes the sphere rotation ω > +30 • and blue ω < −30 • , dark grey is related to the sphere rotation in the range 5 • ≤ ω ≤ 30 • , light grey in the range −30 • ≤ ω ≤ −5 • , medium grey in the range −5 • ≤ ω ≤ 5 • , positive sign means clockwise rotation). Such a colour convention made shear zones clearly observable (only particles within shear zones significantly rotated). There existed a clear grain separation between a clockwise (red) and an anti-clockwise (blue) rotation-the majority of 'red grains' was located within the dominant curve shear zone while the majority of 'blue grains' was placed within the straight radial shear zone (however there existed also a small amount of blue grains within the 'red shear zone' and vice versa). Based on grain rotations, a curved shear zone started to develop along the specimen bottom for the normalized wall translation of u h / h w = 0.02 (Fig. 2Bb). It was fully developed for u h / h w = 0.06. Its thickness was t s = 20 mm (20 × d 50 ). The radial shear zone started later to form for u h / h w = 0.04 (Fig. 2Bd) and for u h / h w = 0.06 connected the curved shear zone. Its thickness was t s = 10 mm (10 × d 50 ). There was a satisfactory agreement between DEM simulation results and real experimental outcomes [41,42] and FE results [19,38].

Helmholtz-Hodge decomposition
The HHD of vector fields is one of the fundamental theorems in fluid dynamics [20,21,43]. It describes a vector displacement increment field in terms of its curl-free and divergence-free components based on potential functions. The unique HHD of the smooth 3D vector field ξ provides the following formula where + ∂ ∂z denotes the divergence operator, ∇× is the curl operator, u denotes the scalar potential field, ν 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 ∇ × ν is called the divergence-free component and is related to vorticity and pure shear (because is incompressible). The harmonic component is related to pure translation.
A variational calculus approach was used [44] which allowed for finding the vector fields ∇u and ∇ × ν by examining the difference between the unknown vector field and provided field ξ (see Eqs. 11,12). By requesting that this difference is minimum (the minimum was found by requesting that the derivatives of the functionals were equal to zero, Eqs. 13,14), the vector fields ∇u and ∇ × ν were explicitly determined. The explicit calculation for ∇u was given in Eqs. 15 and 16 and for ∇ × ν in Eq. 17.
The HHD on irregular grids has been already used e.g. in graphics [43]. In order to create a grid, the centre of each sphere was a node in a Delaunay triangulation [45] and the i-th node had the coordinate r i . Then the discrete piecewiseconstant vector field ξ ( r i ) = k ψ k ( r ) ξ k was created by assigning the constant vector value ξ k to each k-th tetrahedron (ψ k is the piecewise-constant 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 ν are the piecewise linear functions described using a piecewise-linear basis shape function φ i ( r ), their derivatives ∇ will be piecewise-constant, hence the solution for the piecewise-constant ξ( r ) discrete vector field is exact [43].
The decomposition of the discrete vector field ξ was separately solved for each component of u and ν by finding the minimum of 2 following quadratic functionals [43] using the variational calculus principle [44] 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 i (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( ν) is found by requiring that their functional derivatives are zero Solving for u( r ) using the variational calculus : Similarly for ν The integrals in Eqs. 16 and 17 may be 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 is solved for the unknows u i and ν i using standard methods (e.g. inverting matrix or conjugate gradient).
where |T k |-the tetrahedron volume (triangle area), (∇φ i ) k -the 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 ) N (i)-the set of all tetrahedrons (triangles) containing the i-th node.
Equations 18 and 19 describe the i-th row of 2 sparse matrices and were numerically solved for the unknowns u i and ν i using the Eigen library with a bi-conjugate gradient stabilized solver [45]. The third component of Eq. 10 the harmonic vector field (which contains a non-integrable field component) is determined as h = ξ − ∇u − ∇ × ν. The sphere displacement increments were calculated during 10,000 iterations ( u h / h w = 0.002) [16].
In order to obtain a unique solution, appropriate boundary conditions have to be assumed [20,43]. The system of linear equations in HHD was solved using the following general boundary conditions: ∇ × ν (divergence-free component) was tangential to the domain boundary v| ∂T = 0 and ∇u (curl-free 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 can be found in Denaro [46]. Note that a change of these boundary conditions suggested in Petroneto et al. [47] may create an invalid or ill-posed problem [48]. The so-called Hodge-Morrey-Friedrichs boundary conditions may be also used [21]. The boundary conditions obviously influence vector fields close to specimen boundaries. In particular when the spheres' number is low; vortex-structures may be solely detected in the specimen centre since the vector field ∇ × ν is forced by boundary conditions to be parallel to boundaries. In our calculations 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 based on preliminary calculations (e.g. in the case of 20-times smaller number of spheres with the diameter of 10 mm along the specimen height, a very low number of vortices was detected due to a strong effect of boundary conditions). The effect of boundary conditions may be weakened by introducing virtual particles outside the boundary [47]. The calculated vector fields were plotted using the line integral convolution method [49]. No additional smoothing techniques [43] were used in the detection's method to find sinks, sources and vortices.

Passive wall translation
The calculation results are described in Figs. 3, 4, 5 and 6. Figure 3 shows the evolution of the vector displacement increment field ξ during the normalized horizontal wall translation u h / h w (sphere displacement increment directions are marked by white arrows). The scale attached denotes the sphere displacement increment vector length during 10,000 iterations in [mm/iteration] which changes between 0 and 1 mm. Based on the displacement increment vector length and vector direction changes, a curved shear zone between the wall bottom and free upper boundary already started in the first calculation step. Later it moved to the right to reach its ultimate position due to wall friction along the bottom that was close to the maximum resultant normalized horizontal earth pressure force of Fig. 2A (u h / h h = 0.04). Behind the curved shear zone the material was totally rigid. A radial shear zone was approximately created for u h / h w = 0.065.
The evolution of the scalar field gradient ∇u (curl-free component related to compressibility) during normalized wall translation u h / h w is described in Fig. 4. The scale attached denotes the scalar potential u in [mm 2 /iteration]  (sign (−)-dilatancy, sign (+)-contractancy), changing between −4 and 2 mm 2 /iter. The green circles describe the sources (local minima of the scalar potential u -local dilatancy minima) and the red circles denote the sinks (local maxima of the scalar potential u -local contractancy maxima). The local extrema of the scalar field u were defined by requesting that the values of u in all neighbouring nodes in the mesh created by the Delaunay triangulation were smaller/larger than u i at the node in question.
Initially global contractancy and later global dilatancy occurred in the granular specimen. The global dilatancy was the largest after the stress peak (Fig. 4d) and later dimin-ished. The local contractancy maxima and local dilatancy minima started to develop in two main shear zones before the stress peak (Fig. 4c). In the residual state (u h / h w ≥ 0.075, Fig. 2A), local regions of dilatancy and contractancy alternately happened along both shear zones (with the prevalence  h w : a 0025, b 0.025, c 0.04,  d 0.05, e 0.075, f 0.1, g 0.125 and h 0.15 (scale denotes increment vector length h) of dilatancy). This outcome is in accordance with the distribution of local void ratio in a shear zones in calculations by DEM [16]. Note that local dilatant regions may be connected to the collapse of main force chains and creation of vortices and local contractant regions may be connected to the build-up of main force chains and disappearance of vortices [16]. Figure 5 presents the evolution of the vector field curl ∇ × ν (divergence-free component related to vorticity) during normalized wall translation u h / h w . The scale denotes the component of the vector potential ν perpendicular to the specimen in [mm 2 /iteration]), changing from −8 mm 2 /iter up to 8 mm 2 /iter. The green circles describe the local minima of the scalar field v (right-handed vortices) and red circles the local maxima of the scalar field v (left handed vortices). The local extrema of the scalar field v were interpreted as vortices and were calculated in the same way as the local extrema of the scalar field u. The vortex-structures appeared from the wall translation beginning. They were immediately concentrated in regions of the shear zones' occurrence. Thus the ultimate shear zone pattern turned out to be encoded in the grain kinematics from the deformation beginning. This outcome is in accordance with our earlier calculation results for plane strain compression based on displacement fluctuations [18] using the method described in Gould et al. [28] and calculation results based on bottlenecks in force transmission through the contact network [50]. Right-handed and left-handed vortices alternately occurred. The right-handed vortices dominated in the curved shear zone and left-handed vortices dominated in the radial shear zone. Their distance along the shear zones was different.
Finally Fig. 6 shows the evolution of the harmonic increment vector field h during normalized wall translation u h / h w . The scale denotes the increment vector length during 10'000 iterations in [mm/iteration] changing from 0 mm up to 0.8 mm. The evolution h was practically the same independently of u h / h w since the wall continuously moved at the same displacement increment. Figure 7 shows the number of vortex-structures N detected in the entire granular specimen (Fig. 7A), curved shear zone (Fig. 7B) and radial shear zone (Fig. 7C). They were continuously created. Their number rather increased from the beginning up to the residual state (u h / h w = 0.06, Fig. 2A). The number of right-handed vortices was significantly larger in the entire specimen (Fig. 7A). The right-handed vortices were dominant in the curved shear zone (Fig. 7B) and lefthanded ones were dominant in the radial shear zone (Fig. 7C). The maximum/minimum number of right-handed vortices in the curved shear zone was 42/22 (u h / h w > 0.06). The maximum/minimum number of left-handed vortices in the radial shear zone was 13/1 (u h / h w > 0.06).
During the entire wall translation (in the range of u h / h w = 0−0.15), the predominant period of the right-handed vortices was equal to 4% of u h / h w (entire specimen and curved shear zone) (Figs. 8,9). The predominant period of the left-handed vortices was also equal to 4% of u h / h w (entire specimen and radial shear zone) (Figs. 8,9). For the case of Fig. 9Ca, the predominant periods could not be precisely determined. For the residual state (u h / h w = 0.06 − 0.15), the predominant period of the right-handed vortices in the curved shear zone of the left-handed vortices in the radial shear zone was smaller (1-2% of u h / h w , Fig. 10).
The number of sources (local dilatancy minima) and sinks (local contractancy maxima). was larger than of vortices by the factor 3-7 (Fig. 11A). However their predominant period (4% of u h / h w ) was similar (Fig. 11B).

Plane strain compression
Our numerical outcomes with respect to vortex-structures were next checked during quasi-static plane strain compression. The details of DEM 3D calculations are given in Kozicki and Tejchman [18]. The granular specimen used in DEM had the same size as in the experiments by Vardoulakis [51], namely: the width b = 4 cm, height h = 14 cm and depth l = 8 cm (out-of-plane direction) (Fig. 12a). The linear grain distribution curve was assumed; the grain diameter range was between 1.25 and 3.75 mm with d 50 = 2.5 mm. About 56,000 spheres were used with the same material constants. The initial void ratio was e o = 0.53. The flexible vertical walls were assumed to model the membrane surrounding the specimen in experiments (Fig. 12a). Both the front and rear specimen sides 4 cm 2 × 14 cm 2 were blocked in a perpendicular direction to the specimen to enforce plane strain conditions. The bottom surface 4 cm 2 × 8 cm 2 was fixed in a vertical direction and the top surface 4 cm 2 × 8 cm 2 was subjected to a constant vertical displacement u ν . 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. The evolution of the mobilized internal friction angle (calculated with principal stresses from the Mohr's equation) versus the vertical normal strain and displacement vectors of spheres in the specimen are shown in Fig. 12b, c. Similarly as in real experiments [51], the initially dense specimen showed an asymptotic behaviour; it exhibited initially small elasticity, hardening (connected first to contractancy and then dilatancy), reached a peak strength at about of ε 1 = 5%, gradually softened and dilated reaching a residual state at the large vertical strain of 25-30% (Fig. 12b). During deformation a distinct internal inclined shear zone occurred inside the sand specimen which was marked by shear strain, larger grain rotation and volume increase. The thickness of the inclined interior shear zone t s was on average in the residual state for ε 1 = 30% about t s = 25 mm (10 × d 50 ) based on strain deformation in the specimen. The calculated shear zone inclination to the bottom was 60 • at ε 1 = 10% and 67 • at ε 1 = 30%. Based on the cumulative grain rotation, the internal inclined shear zone might be noticed for ε 1 ≈ 3% [18]. Due to a rather small number of spheres along the specimen width (about 16), the effect of boundary conditions during calculations of vortex-structures was weakened by introducing virtual particles outside boundaries [47] (Fig. 12c). Artificial nodes were added in the Delaunay's triangular mesh at the distance of up to 50 mm around the specimen (with the grid inter-node distance of 1 mm) (Fig. 12c). The vector ξ in these artificial     1 = 8% (marked in red) [18], b mobilized internal friction angle φ versus ε 1 [18] and c displacement vectors of specimen spheres and artificial nodes around specimen for ε 1 = 10% (vectors are enlarged by factor 200) (colour figure online) nodes was calculated using the Gaussian averaging for true specimen nodes with the averaging radius of 80 mm (2 × b). Some wild vectors in Fig. 12c appeared if single grains suddenly undergo large displacements. Figure 13 presents the evolution of the vector field curl ∇× ν (divergence-free component related to vorticity) during the normalized top displacement u ν / h in the vertical mid-section slice with the area of 4 cm 2 × 14 cm 2 and thickness of 5 × d 50 (1.25 cm) (Fig. 12a) cut out from the granular specimen Fig. 13 Evolution of vector field curl ∇ × ν (divergence-free component related to vorticity) 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% [scale denotes component of vector potential ν perpendicular to specimen in (mm 2 /iteration)], green circles describe local minima (right-handed vortices) and red circles local maxima (left handed vortices) (colour figure online) [18]. The vortex-structures appeared from the deformation beginning. They were immediately concentrated in the region of the mean shear zone occurrence. Right-handed vortices were dominant in the shear zone.
In order to fully validate the proposed method for detecting vortex-structures, several various boundary value problems with shear localization will be investigated as. e.g. during simple and direct shearing and wall shearing [38]. The grain properties, shear rates and grain size distribution ranges will be changed. Moreover, 3D vortex-structures will be calculated during plane strain compression [18] using HHD. The effect of the different definition of the rolling velocity [33] will be also investigated. In parallel, similar calculations will be carried out within micro-polar hypoplasticity [38].

Conclusions
The following conclusions may be listed from DEM simulations of 2D quasi-static granular vortex-structures in sand the mean grain diameter of 1 mm during a quasi-static passive earth pressure problem: -The HHD allowed for separating a vector field into the sum of three uniquely defined components: curl free, divergence free and harmonic. It proved to be an objective, universal and effective technique for identifying all vortex-structures during granular flow which was directly based on single grain displacement increments (but not on displacement fluctuations). The method did not use any additional non-objective parameters. A large number of spheres was however required to avoid the effect of boundary conditions assumed. The size of vortexstructures could not be deduced since they corresponded to points only which were associated with the centre of shear zones. The method may be used for the detection of 3D vortices. -A strong connection between the location of vortexstructures and progressive shear localization was found out. The vortex-structures were the precursor of shear localization since they clearly concentrated in the area where a curved and radial shear zone ultimately later formed. Thus the ultimate shear zone pattern was 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 a reliable indicator of shear localization. They developed from the deformation process beginning. They solely emerged in main shear zones. They had a tendency to move along shear zones. Their number varied and 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 curved shear zone, the predom-inant period of right-handed vortices was 4% of u/ h during the entire wall movement. In the radial shear zone, the predominant period of left-handed vortices was also 4% of u/ h. -In the residual state, local regions of dilatacy and contractancy alternately happened along globally dilatant shear zones with a dominance of local dilatancy. -An early prediction possibility of shear localization through vortex-structures may open new perspectives for a detection of impending failure in granular bodies (inherently connected with shear localization) within continuum mechanics.