Investigations of formation of quasi-static vortex-structures in granular bodies using DEM

The paper presents some two-dimensional simulation results of vortex-structures in cohesionless initially dense sand during 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 flow displacement field from DEM calculations was used. This approach enabled us to distinguish both incompressibility and vorticity in the granular displacement field.


Introduction
Granular vortex-structures (swirling motion of several grains around its central point) were frequently observed in experiments on granular materials [1], [2] and in calculations using the discrete element method (DEM) [3]- [6]. They became apparent when the motion associated with uniform (affine) strain was subtracted from the actual granular deformation. A dominant mechanism responsible for the vortex formation was buckling of force chains [3], [5]. The collapse of main force chains lead to a formation of larger voids and their build-up to the formation of smaller voids [3]. The vortex-structures have been mainly observed in shear zones [4] which are the fundamental phenomenon in granular bodies [7].
The objective of this paper is to report the results of comprehensive 2D studies by DEM on quasi-static vortex-structures in sand behind a rigid wall during its passive translation by using the Helmholtz-Hodge decomposition of a displacement vector field [8], [9]. This detection method has the following advantages over other techniques [4]: 1) it operates directly on displacement vectors, 2) it is designed for both 2D and 3D flow and 3) it does not need any additional parameters for calculations. DEM analyses were carried out with spheres with contact moments [3], [4], [10] to approximately capture the irregular grain shape (the contact moment increments were calculated by means of the rolling stiffness multiplied by the angular rotational increment vectors). In order to accelerate the computation time, some simplifications were assumed in numerical studies: larger spheres, linear grain size distribution and linear normal contact. Attention was paid to a relationship between the vortex occurrence and shear localization. In addition, the predominant periods of vortex-structures were calculated.

DEM results of passive earth pressure model tests
In order to simulate the behaviour of real sand, the 3D spherical discrete model YADE, developed at University of Grenoble [11] was used by taking advantage of the socalled soft-particle approach (i.e. the model allows for particle deformation which is modelled as their overlap).
The DEM calculation results were described in detail in [3]. The quasi-static simulations were performed for a 2D initially dense sand body of 0.40 m length and 0.20 m height in order to compare with experiments with Karlsruhe sand (mean grain diameter d 50 =0.5 mm) [7]. The vertical retaining wall and the bottom of the granular specimen were assumed to be stiff and very rough. Since the experiments were idealized as a 2D boundary value problem, 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 particles with d 50 =1.0 mm, characterized by a linear grain size distribution, were assumed (grain size range was 0.5-1.5 mm, 62'600 particles). The initial void ratio of sand was e o =0.62 (volumetric weight J=25.5 kN/m 3 ). Figure 1A shows the evolution of the resultant normalized horizontal earth pressure force (earth pressure coefficient) K p =2E h /(γh 2 d 50 ) versus the normalized horizontal wall displacement u/h (h=0.2 m, E h -the horizontal force acting on the wall) for d 50 =1 mm. The specimen exhibited the initial strain hardening up to the peak (u/h=0.038), followed by some softening before a common asymptote was reached. The earth pressure coefficient was K p max =30 for d 50 =1 mm. The earth pressure coefficient was K p max =30 for d 50 =1 mm. It can be thus anticipated that for d 50 =0.5 mm (real sand), K p max should be about 25-27. The value of K p max =30 for d 50 =1 mm was slightly smaller than K p max =31 obtained by FEM (d 50 =0.5 mm) [7] and was closer to engineering earth pressure coefficients [7].
The distribution of single sphere rotations Z during wall translation is presented in Fig.1B (red denotes the sphere rotation Z>+30 o and blue Z<-30 o , dark grey is related to the sphere rotation in the range 5 o dZd30 o and light grey to the range -30 o dZd-5 o , positive sign means clockwise rotation). Accepting such a colour convention makes shear zones clearly observable (only particles within shear zones significantly rotate). There exists a clear grain separation regarding clockwise (red) and anti-clockwise (blue) rotation. Except of two main shear zones (curved and radial shear zone), some secondary shear zone might also form [3]. The curved shear zone started to develop along the specimen bottom for the normalized wall translation u/h=0.02 (Fig.2Bb). It was fully developed for u/h=0.06. Its thickness was t s =20 mm (20×d 50 ). The radial shear zone began to form for u/h=0.04 (Fig.1Bd) and for u/h=0.06 connected the curved shear zone. Its thickness was t s =10 mm (10×d 50 ). The horizontal force strongly fluctuated at the residual state [12]. There was a good qualitative agreement between DEM simulation results and real corresponding experimental outcomes [7], [13].

Numerical results using Helmholtz-Hodge decomposition
The Helmholtz-Hodge decomposition (HHD) of 3D vector fields is one of the fundamental theorems in fluid dynamics [8], [9]. It describes a vector displacement field in terms of its curl-free and divergence-free components based on potential functions: , wherein is the gradient, denotes the divergence operator, is the curl operator, u denotes the scalar potential field, is the vector potential field and denotes the harmonic vector field. The gradient of the scalar potential function 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 (because is incompressible). The harmonic component is related to pure translation. A discrete vector field decomposition was done according to the method described in [8], [9].
The system of linear equations was solved using the following general boundary conditions: was tangential to the domain boundary =0 and was orthogonal to the boundary domain =0. The calculations were carried out for one granular specimen. Figure 2 presents the evolution of the vector field curl (divergence-free component related to vorticity) during normalized wall translation u/h. The green circles describe the local minima (right-handed vortices) and the red circles the local minima (left handed vortices). Vortex-structures appeared from the wall translation beginning. They were immediately concentrated in shear zones only. Right-handed and lefthanded vortices alternately occurred. The right-handed vortices dominated in the curved shear zone and lefthanded vortices dominated in the radial shear zone. Their distance was different. The maximum/minimum number of left-handed vortices in the radial shear zone was 24/1. During the entire wall translation (u/h=0-0.15), the predominant period of the right-handed vortices was equal to 1.5% of u/h (entire specimen) and to 1.5% of u/h (curved shear zone). The predominant period of the left-handed vortices was equal to 2% of u/h (entire specimen) and to 4% of u/h (radial shear zone). For the residual state (u/h=0.06-0.15), the predominant period of the righthanded vortices in the curved shear zone was 1.25% of u/h and of the left-handed vortices in the radial shear zone was 1% of u/h. Figure 3 shows the vector displacement field during the normalized horizontal wall translation u/h. Sphere displacement directions are marked by white arrows. Based on the displacement vector length and direction changes, a curved shear zone between the wall bottom and free upper boundary already began in the first calculation step. Later it moved to the right to reach its ultimate position (Fig.1Bd)  The evolution of the scalar field gradient (curlfree component related to compressibility) during normalized wall translation u/h is described in Fig.4. The green circles describe the sources (local dilatancy minima) and red circles the sinks (local contractancy maxima). Initially global contractancy and later global dilatancy occurred in the granular specimen. The dilatancy was maximum after the peak and later it diminished. In the residual state (u/ht0.075, Fig.1A), local regions of dilatancy and contractancy alternately happened along both shear zones (with the prevalence of dilatancy) (Fig.4).

Conclusions
The following main conclusions may be listed from DEM simulations of vortex-structures in sand during a passive earth pressure problem: A strong relationship was found between the location of vortex-structures and shear localization. Vortexstructures were the precursor of shear localization since they clearly appeared at the place of a curved and radial shear zone from the beginning of the deformation process, i.e. significantly earlier than e.g. based on single grain rotations. They vortex-structures solely appeared in shear zones during the entire deformation process. 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.
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 predominant period of right-handed vortices was 1.5% of u/h during the entire wall movement (calculated based on the number of vortices versus u/h using FFT). In the radial shear zone, the predominant period of left-handed vortices was 4% of u/h. In the case of the residual state (u/h=0.06-0.15), the predominant period of right-handed vortices in the curved shear zone was 1.25% of u/h and of left-handed vortices in the radial shear zone was 1% of u/h.
In the residual state, local regions of dilatacy and contractancy alternately happened along main shear zones with a dominance of local dilatancy.