3D DEM simulations of monotonic interface behaviour between cohesionless sand and rigid wall of different roughness

The paper deals with three-dimensional simulations of a monotonic quasi-static interface behaviour between cohesionless sand and a rigid wall of different roughness during wall friction tests in a parallelly guided direct shear test under constant normal stress. Numerical modelling was carried out by the discrete element method (DEM) using spheres with contact moments to approximately capture a non-uniform particle shape. The varying wall surface topography was simulated by a regular mesh of triangular grooves (asperities) along the wall with a different height, distance and inclination. The calculations were carried out with different initial void ratios of sand and vertical normal stress. The focus was to quantify the effect of wall roughness on the evolution of mobilized wall friction and shear localization, also to specify the ratios between slip and rotation and between shear stress/force and couple stress/moment in the sand at the wall. DEM simulations were generally in good agreement with reported experimental results for similar interface roughness. The findings presented in this paper offer a new perspective on the understanding of the wall friction phenomenon in granular bodies.


Introduction
Soil-structure interfaces are frequently encountered in geotechnical engineering, e.g. foundations, tunnels, retaining walls, anchors, silos, piles and geotextiles. They play a major role in the interaction between soils and structures with respect to a static, dynamic and fatigue mechanical behaviour and durability performance. Interface mechanical properties are affected by properties of both the contacting soil and opposing interface. Therefore, a robust understanding of soil-interface behaviour is essential for geotechnical designs. The interface behaviour is characterized by the formation of a wall shear zone with a certain thickness in the soil adjacent to the structure, i.e. a thin zone of intense shearing with both pronounced grain rotations and volume changes. The determination of the thickness of the wall shear zone is of major importance for estimating the shear resistance and forces transferred from the surrounding soil to the structure, resulting in evaluating of the structure strength. The thickness of a wall shear zone depends on several factors such as wall roughness and stiffness, initial density and mean grain diameter of soil and shearing velocity [49,53,55]. In problems involving the interface behaviour, Coulomb's friction law is usually used, based on the assumption of a constant ratio between the shear and normal stresses on the interface. However, such assumption provides solely an approximate description of the soil-structure interface since: (1) the mobilized friction angle can change significantly during shearing along rough or very rough walls in contact with initially dense granulate [49,53] and (2) the mobilized wall friction angle is not a state variable, as it depends on a number of factors, such as boundary conditions, contact pressure level, initial stress state and specimen size [49,53]. Therefore, the laboratory wall friction angles obtained in a test configuration (in particular, the peak values) cannot always be directly transferred to other boundary value problems [49].
The objective of the current paper is to carefully study a monotonic quasi-static interface behaviour between cohesionless sand and a rigid wall of different surface topography in wall friction tests using a direct shear box under conditions of constant normal stress. The different artificial wall surface topography was created by regularly arranged triangular grooves (asperities) at the same spacing in the form of a standard saw-tooth surface. The simulations were carried out with the discrete element method (DEM) under 3D conditions. The particle-based open-source code YADE, elaborated at the University of Grenoble, was used for DEM simulations. The effects of the height, distance and inclination of grooves (asperities) were carefully studied for the varying initial void ratio of sand and pressure. Some calculation results were directly compared with our corresponding monotonic quasi-static wall friction tests on cohesionless 'Karlsruhe' sand in a parallelly guided direct shear box, performed at the Karlsruhe University [49,53]. To approximately simulate the irregularity of sand particles' shape of 'Karlsruhe' sand, spheres with contact moments were used. The focus was on the effect of wall roughness related to two aspects: (1) the evolution of both mobilized wall friction and wall shear zone and (2) the distribution of grain displacements and rotations, shear stresses/forces and couple stresses/moments in the granular assembly directly at the wall.
The paper includes a few novel points: (1) comprehensive analyses of the interface behaviour with the real mean grain diameter of sand by taking into account the effect of different wall roughness parameters (height, distance and inclination of grooves), initial void ratio of sand and pressure effect on the wall friction characteristics, (2) 3D DEM simulations of the particle assembly (it is well known that there exist some fundamental discrepancies between numerical predictions for models of 2D and 3D particles), (3) the determination of displacements and rotations, shear stresses and couple stresses, wall forces and wall moments in the granular specimen at the walls of different roughness and (4) the proposal of wall boundary conditions for micropolar continua. The limitations of the current DEM study are related to two issues: (1) the approximate shape of 'Karlsruhe' sand grains was assumed and (2) the experimental wall roughness was not faithfully reproduced. The findings presented in this paper can help to better understand the wall friction mechanism, to design geotechnical systems with enhanced strength and to define wall boundary conditions within models of micro-polar continua.
The paper is arranged as follows. A brief summary of the related past work was given in Sect. 2. In Sect. 3, the proposed numerical DEM framework was presented in detail. The model calibration was discussed in Sect. 4. Results of the effect of wall roughness on wall friction angle and shear localization with key findings were described in Sects. 5. The main results were summarized and some conclusions were stated in Sect. 6.

Literature overview
The interface between granular material and structure has been investigated using various testing devices and methods [53], e.g. direct shear apparatus [2,7,8,17,19,40,42,43,46,48,49,64,71,73], torsional ring shear apparatus [21,33,37,69], ring shear device [3,23], simple shear apparatus [57,58], plane strain apparatus [49], Couette apparatus [1,33], wear tester [16], three-dimensional simple shear apparatus [11], ring simple shear apparatus [31] and in experiments with piles [61], anchors [65] and silos [49]. The experimental results showed a pronounced effect of the wall roughness, grain size, grain distribution, pressure level, initial density, specimen size and velocity on the peak wall friction angle and wall shear zone thickness. The shear zone thickness was found to increase with increasing wall roughness, grain size, pressure, shear strain rate, specimen size and to decrease with increasing initial unit weight. The mobilized wall friction angle at the peak grew with increasing wall roughness, grain size, initial unit weight and velocity, and decreased with growing pressure and specimen size. Moreover, large void fluctuations, grain mixing and grain segregation were observed in the wall shear zone. The maximum interface strength was achieved, e.g. for an asperity distance to mean grain diameter ratio between 1.0 and 3.0, and an asperity height to mean grain diameter ratio greater than 0.9 [8]. An asperity angle of 50 o or greater yielded the maximum efficiency for any given asperity spacing or height [8]. The maximum wall friction angle was found to be larger than the internal friction angle of soils due to the passive resistance caused by surface asperities [13,49]. It always varied in a bi-linear fashion as a function of the normalized roughness [19,40,48]. The critical normalized roughness was mainly about 0.4-0.5, and above this value, the maximum wall friction angle insignificantly increased. The pure slip during wall shearing was registered by Uesugi [57,58], increasing with decreasing wall roughness. In these experiments, the ratio between wall grain rotations and their slips was not measured.
Besides the experimental studies, several numerical DEM and FEM analyses were carried out to investigate the interface behaviour in granular materials. Modelling the interface thickness within continuum mechanics using FEM can be only performed with constitutive models possessing a characteristic length of microstructure [9,10,20,49,50,53,55,59,60]. Moreover, the constitutive models have to take the salient behaviour of granular materials into account by considering major influential factors such as the initial density, pressure sensitivity and mean grain diameter of soils [52,54]. The boundary conditions at interfaces with consideration of a characteristic length of microstructure were investigated following different enriched approaches, e.g. within micropolar elasto-plasticity [10,49,53,59], strain gradient elasto-plasticity [60] and micro-polar hypoplasticity [9,20,50]. Most of the calculations were carried out with very rough interfaces. Different micro-polar boundary conditions were proposed in [9,50,55] for describing the wall roughness. In [55], the boundary conditions along the horizontal rigid wall suggested inclusion of two ratios connected to the normalized wall roughness (a ratio of the micro-polar rotation multiplied by the mean grain diameter and the horizontal displacement and a ratio of the horizontal shear stress multiplied by the mean grain diameter and the horizontal couple stress). To better understand microscopic phenomena during wall friction, DEM calculations were also carried out [4,5,12,14,15,22,24,25,62,63,70,77,78]. There exist many numerical studies of wall friction using DEM under 2D conditions (e.g. [14,15,22,62,63,77,78]) and only a few under 3D conditions [4,5,12,24,25,70]. To facilitate the interpretation of macroscale responses, microscale metrics such as contact normal force distribution, contact networks, mobilization of friction, and particle rotation were calculated to elucidate the wall friction mechanism. The DEM simulations also exhibited a bi-linear relation between the interface resistance and normalized wall roughness [25,70]. The critical normalized roughness was found to be about 0.4 [25,77]. However, no effort is known to us that was performed in experiments and numerical DEM simulations on wall friction to determine the ratios between grain rotations and slips, shear stresses and couple stress and forces and moments at the wall of different roughness. These ratios are of importance for defining wall boundary conditions within models of micro-polar continua [49,53,55].

3D DEM model
To evaluate the effects of soil properties and interface roughness on the behaviour of real sand, the 3D spherical discrete element model YADE developed at the University of Grenoble [27,45,68] was employed. To approximately simulate the irregularity of the particles' shape of 'Karlsruhe' sand, spheres with contact moments were assumed [29,30,39,66]. DEM has natural predisposition to account for the material non-uniformity as complex global constitutive relationships are replaced by simple local contact laws in DEM. 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 strain localization. DEM may be used to frictional [29,30,39,66] and frictional-cohesive materials [32,38,72]. The disadvantages are related to an enormous computational cost and an extensive calibration based on experimentally measured macro-scale properties. The algorithm used in the present DEM is based on a description of particle interactions in terms of force laws and involves in general two main steps. First, interaction forces between discrete elements are computed, based on constitutive laws. Second, Newton's second law is applied to determine for each discrete 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 (interpreted as a local contact deformation). The role of the particle shape was highlighted in [28,74,75]. A linear elastic normal contact model was used only. In compression, the normal force was not restricted and could increase indefinitely. Figure 1 shows the mechanical response of the contact model when using spheres with contact moments. The DEM model can be summarized as follows [27,45,68]: jjF s jj À jjF n jj Â tan l 0; ð4Þ where F n -the normal contact force, U-the overlap between discrete elements, F s -the tangential contact force, F s;prev -the tangential contact force from the previous iteration, Ñ-the unit normal vector at each contact point, X s -the relative tangential displacement of the sphere centre, K n -the normal contact stiffness, K s -the tangential contact stiffness, l-the Coulomb inter-particle friction angle, R-the element radius, R A and R B -the contacting grain radii, E c -the elastic modulus of the grain contact, m c -the Poisson's ratio of the grain contact, M-the contact moment, K r -the rolling contact stiffness, bthe dimensionless rolling stiffness coefficient, x-the resultant angular rotation between two elements, g-the dimensionless limit rolling coefficient, F k damp and M k dampthe damped contact force and moment, F k and M k -the kth components of the residual contact force and contact moment vector, m k and _ x ! k are the kth components of the translational and rotational velocities of spheres and a dthe positive numerical damping coefficient smaller than 1 [6] (sgn(•) returns the sign of the kth component of the translational and rotational velocity). 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 m [28,30]. The effect of damping was negligible in quasi-static calculations [28,30].
The five main local material parameters are necessary for our DEM simulations: E c (modulus of elasticity of the grain contact), m c (Poisson's ratio of the grain contact), l (inter-particle friction angle), b (rolling stiffness coefficient) and g (limit rolling coefficient). In addition, a particle radius R, particle mass density q and numerical damping parameter a are required. The DEM material parameters: E c , m c , l, b, g and a were calibrated using the corresponding homogeneous axisymmetric triaxial laboratory test results on Karlsruhe sand with the different initial void ratio and lateral pressure [26,67]. The procedure for determining the material parameters in DEM was described  [15,24,62] in detail by Kozicki et al. [28,30]. Note that the representative elastic contact moduli E c and m c are different from the elastic moduli of grains.
The DEM results were directly compared with the corresponding laboratory tests on 'Karlsruhe' sand in a parallelly guided shear device under constant vertical normal stress [49,53]. The sand specimen size was 100 9 100 9 20 mm 3 . The sand was initially dense or initially loose. The vertical normal stress was varied between 50 and 200 kPa. The experiments were carried out with three different wall roughnesses R max . We classified the wall roughness as smooth (0 \ R max B 0.1 9 d 50 ), rough (0.1 9 d 50 \ R max \ 0.5 9 d 50 ) and very rough (R max C d 50 ), where R max -the maximum vertical distance between peaks and valleys evaluated along the wall over the length of 3 9 d 50 [55,57]. The rough wall was obtained in a corrosion chamber and a very rough wall was obtained with the aid of random glueing of particles of 'Karlsruhe' sand with the mean diameter equal to and higher than 0.5 mm to the wall. The index properties of 'Karlsruhe' sand are: the mean grain diameter d 50 = 0.50 mm, grain size between 0.08 mm and 1.8 mm,  [27] for different initial lateral pressures: r c = 50 kPa, r c = 200 kPa and r c = 500 kPa Table 1 Material parameters assumed in all DEM simulations

Material parameters Value
Modulus of elasticity of grain contact E c (MPa) 300 Poisson's ratio of grain contact v c (-) 0.3 Inter-particle friction angle l (8) 1 8 void ratio e max = 0.84. The sand grains were classified as sub-rounded/sub-angular.

Model calibration based on triaxial compression
A triaxial compression test is the most frequently geotechnical test used for calibration of soils. To determine the material parameters in DEM (E c , m c , l, b, g and a), a series of numerical homogeneous quasi-static triaxial compression tests with rigid smooth walls on cohesionless sand were initially performed [28,30]. The DEM results were compared with corresponding comprehensive experimental triaxial compression results with 'Karlsruhe' sand for the different initial void ratios and lateral pressures [26,67]. For simulations, a cubical specimen of size 10 9 10 9 10 cm 3 , composed of about 8000 spherical particles with contact moments was constructed. The grain diameter of sand linearly varied between 2.5 mm and 7.5 mm and its mean grain diameter was d 50 = 5 mm (10 times larger than the real one). The mass density was 2600 kg/m 3 . The granular assembly was prepared by putting spheres of a random radius according to the grain distribution curve (without gravity) into a cubical container with six external walls, which had a regular cubical grid with a particle distance of 10 mm. In order to obtain a desired initial density owing to grain overlapping, the interparticle friction angle was varied between 0°and l (initially dense sand) and between 89°and l (initially loose sand) to exactly reproduce the target initial void ratio. During dynamic compression to the desired confining pressure r c , grains bounced against each other and moved in random directions; thus, their initial ordered arrangement became entirely random. The assembly was then allowed to settle to a state where the kinetic energy was negligible and then friction coefficient was set to the target inter-particle friction angle l. The DEM simulations were carried out for initially dense sand (initial void ratio e o = 0.53) for three different lateral pressures r c = 50 kPa, 200 kPa and 500 kPa ( Fig. 2B and C) with the parameters listed in Table 1. A satisfactory agreement was obtained, in particular, for the stress-strain curves (Fig. 2B). A comprehensive comparison between the DEM and experimental results was discussed in [28,30]. The comparative DEM calculations of a triaxial compression test were also performed with a nonlinear contact law [36] following Hertz [18] and Mindlin and Deresiewicz [35]. They demonstrated small differences as compared to the results of linear contact law.   Table 2 Calculated values of peak wall friction angle u w,max , residual wall friction angle u w,res and residual volumetric strain e v for different normalized wall roughness parameter R n (regular grooves)

DEM model and preliminary simulations
Monotonic wall shear tests in a direct shear box were simulated with DEM. The 3D granular specimen included 80,000 spheres with contact moments The advantages of a direct shear box test in measuring soil properties are: the simplicity of both the system, specimen preparation and testing procedure. The shortcomings are: deformation and stress fields are non-uniform within the box, the interface area may diminish during shearing, principal stresses are not known, shear strength is larger than the one from triaxial tests or simple shear tests, stress concentrations occur at ends, initiating slip failure and then the shear-induced displacement cannot be separated from the contact slip [49,51]. The 3D DEM model is shown in Fig. 3. The specimen length (l = 100 mm) and the height (h = 20 mm) ( Fig. 3A) were the same as in the experiment [49,53]. To prevent locking of particles at the bottom corners during shearing, the gap equal to the maximum grain diameter was left between the bottom and vertical walls (see the zoom in Fig. 3A) as in the experiment [49,53]. The sand leakage during shearing had a minor effect on void ratio and volumetric strain (\ 1%). The maximum number of spheres beyond the box after the test was solely 400. The comparative calculations without a gap indicated similar results; however, the evolution of the wall friction angle showed more fluctuations due to grain interlocking at ends.
The grain diameter of sand linearly varied between 0.25 mm and 0.75 mm with the mean grain diameter of d 50 = 0.5 mm (as in the experiment). All walls confining  The topography of the wall roughness is in the reality random and complex. In DEM simulations, the wall roughness was simulated in various way. A regular rough surface was made of overlapped particles with a different centre distance [12] or particles with the same diameter [15]. A regular saw-tooth surface with the same groove inclination was assumed in [25]. The most realistically was the numerical wall roughness described in [63]. Regular saw-tooth surfaces with the varied asperity height, asperity width and spacing between asperities were chosen. In addition, a non-regular saw-tooth surface was analysed [63]. We divided our DEM simulations on wall friction into two steps. In the first step (current paper), the bottom wall of the direct shear box had the artificial surface roughness, created by regularly arranged triangular grooves (asperities) in the form of a regular saw-tooth surface wherein the triangular grooves had the same distance but a different inclination. In the next step, the DEM simulations will be carried out regular triangular grooves of the same inclination but the different spacing as in [25,63]. The wall roughness was characterized by the normalized wall roughness parameter R n = h g /d 50 , where h g is the groove height and d 50 denotes the mean grain diameter [53,57,58]. The parameter R n was 2.0, 1.0, 0.75, 0.50, 0.25, 0.10 and 0.01 (Fig. 3c). The basic groove distance s g was always the same (s g = 29d 50 ). The groove inclination to the bottom a g diminished with decreasing R n (e.g. a g-= 45 o for R n = 1.0) (Fig. 3C). Some DEM simulations were also carried with the different groove inclination to the bottom a g and groove distance s g for R n = 1.0 (Sect. 5.4). The constant uniform vertical pressure r n was applied to the top area of the shear box. The horizontal velocity of the shear box was small enough to consider the test as quasi-static (the inertial number I was kept below 10e -4 ). The wall friction angle between particles and grooves was assumed to be l sw = l = 18 o . Figure 4 shows the effect of 3D calculations as compared to 2D ones (Fig. 4A) and the effect of the different specimen width D = 2.5-100 mm (Fig. 4B) on the evolution of the mobilized wall friction angle u w = arctan(T/  Figure 4A shows that the 3D simulations significantly diminished the particle oscillation in 2D analyses. The full width of the sand specimen of D = 100 mm was not also needed to be considered (Fig. 4B). Therefore, to strongly reduce the computation time, the specimen width D along the coordinate 'Z' (Fig. 3A) was assumed to be equal to D = 5.0 mm (10 9 d 50 ) in all DEM simulations instead of D = 100 mm as in the experiment.

Effect of initial void ratio of sand and vertical pressure
To analyse the effect of initial void ratio e 0 on the specimen behaviour along the rough bottom wall, two different initial void ratios were assumed: e 0 = 0.55 (initially dense sand) and e 0 = 0.80 (initially loose sand) with the normalized wall roughness parameter R n = 1.0 and vertical pressure r n = 100 kPa. The DEM results with e 0 = 0.55 and e 0 = 0.80 were compared with our experiments [49,53] (Fig. 5). The DEM calculations were also carried out with three different vertical pressures r n : r n = 50 kPa, r n = 200 kPa and r n = 500 kPa (R n = 1.0 and e o = 0.55) (Fig. 6). The evolutions of u w = f(u x ) and e v = f(u x ) (Figs. 5 and 6) with R n = 1.0 are typical for the sand behaviour during a direct shear test [49,53]. Initially, the mobilized wall friction angle grew until it reached a peak value for the displacement of about u x = 0.5 mm and next exhibited softening. The calculations yielded the wall friction angle u w,max = 49.68 at the peak (u x = 0.5 mm) and u w,res = 368 (u x [ 3 mm) at the residual state for initially dense sand, and u w,max = 428 at the peak (u x = 1.5 mm) and u w,res-= 348 (u x = 3 mm) at the residual state for initially loose sand (Fig. 5A). The latter indicated little softening due to the initial densification caused by the vertical pressure r n . Globally, the initially dense sand dilated and the initially  (Fig. 5B). The DEM calculation results showed an acceptable agreement with the numerical results (Fig. 5b) by taking into account that the wall roughness and particle shapes in experiments were not exactly reproduced in DEM simulations. The calculated displacements corresponding to u w,max were smaller by the factor 2 than in the experiment. The calculated residual wall friction angles u w,res were also smaller than the experimental values. The calculated volume changes were too large by factor 2 for e o = 0.55 as compared to the experiments. The peak wall friction angle diminished and the corresponding displacement increased with increasing pressure r n (Fig. 6) in agreement with the experiment [49]. The residual (critical) wall friction angle for R n = 1.0 was almost the same, independently of the initial void ratio and normal pressure (Figs. 5 and 6). The outcome with respect to the effect of initial void ratio on u w,res is in agreement with our both wall friction experiments on very rough wall (Fig. 5) and a pure sand shear test [49]. The outcome with respect to the effect of pressure on u w,res matches pure sand shear test results [49]. The residual volumetric strain reduced with growing pressure (Fig. 6B) as in the experiment [48] and in DEM simulations [15] (wherein the critical void ratio was found to decrease linearly with increasing normal stress).

Effect of wall roughness
The influence of the wall roughness on the sand behaviour was analysed in a series of tests with the different normalized wall roughness parameter R n that varied between R n = 0.01 and R n = 2.0 for initially dense sand (e 0 = 0.55) and vertical pressure of r n = 100 kPa (Fig. 7). The values of the peak wall friction angle u w,max , residual wall friction angle u w,res and volumetric strain e v are given in Table 2. The relationship between the computed values of u w,max and u w,res and R n is shown in Fig. 8. The DEM results were compared with our experiments [53] in Fig. 9.

Evolution of mobilized wall friction angle and volumetric strain
The peak u w,max and the residual wall friction angle u w,res rapidly increased with increasing roughness parameter up to a particular value of surface roughness (called the critical surface roughness) as in the laboratory tests by Hu and Pu [19] and Su et al. [48] and DEM simulations by Jing et al. [25] and Zhang and Evans [70]. In our analyses, the critical surface roughness was equal to R n(crit) = 0.50-0.75, beyond which their effect became negligible. For R n-C 0.75, the wall friction behaviour of sand was similar ( Fig. 7A and B, Table 2). The volumetric strain of sand was also the same for R n C 0.75 and reduced with decreasing R n (Fig. 7B). The relationship between the values of u w,max and R n and between u w,res and R n was bilinear as, e.g. in experiments [19,48] and DEM analyses [25,70] (Fig. 8). The value of u w,max = 188 for R n = 0.01 was obviously equal to the assumed wall friction angle between particles and grooves l sw . A satisfactory qualitative agreement with experiments was achieved (Fig. 9). The slightly lower values of R n(crit) were obtained in the experiments by Hu and Pu [19] (R n(crit) = 0.4) and DEM calculations by Jing et al. [25] (R n(crit) = 0.375) and Zhang and Evans [70] (R n(crit) = 0.4). The similar value of R n(crit) was obtained in the tests by Su et al. [48] (between 0.5 and 1.0).

Distribution of particle rotation, particle displacement and void ratio
The effect of the normalized wall roughness parameter R n on the sphere rotations x in the entire sand specimen is shown at the residual state (u x = 7 mm) in Fig. 10. The values of rotations x were calculated from a cubic averaging cell of the size 5d 50  sphere rotations x and void ratio e across the normalized specimen height h/d 50 at the specimen mid-region at the residual state. The values of u x , x and e were calculated from the averaging cell of the size 5d 50 9 5d 50 9 1d 50 (length 9 width 9 height) being moved by d 50 . Figure 13 presents the relationship between the thickness of the wall shear zone and the normalized roughness parameter R n . The effect of the initial void ratio of sand e o and vertical pressure r n on the distribution of sphere rotations x across the normalized specimen height h/d 50 is demonstrated in Fig. 14 at the residual state. The distribution of the ratio between the sphere rotation multiplied by the mean grain diameter xd 50 and the sphere slip u [(xd 50 )/u ] across the specimen height at the specimen mid-point is shown in Fig. 15 at the residual state. The wall grain rotation-wall grain slip ratio (xd 50 )/u against the inversed normalized wall roughness parameter 1/R n at the grooves' height is given in Fig. 16.
During sand-wall shearing along the bottom with the different normalized wall roughness parameter R n , an almost horizontal dilatant shear zone along the bottom was created, based on particle rotations and increasing void ratio, which are the best indicators for shear localization [49,53]. The thickness of the wall shear zones was based on an inflexion point in the distribution of sphere rotations x where the rotation was x B 5% of the maximum sphere rotation x max in the wall shear zone (Fig. 12B). The thickness of the wall shear zone t s increased in an approximate bi-linear way with growing R n up to R n \ 0.75 only (Figs. 10, 12B and 13). It was: t s = 14 9 d 50 (R n-= 0.75-2.0), t s = 10 9 d 50 (R n = 0.50), t s = 69d 50 (R n-= 0.25), t s = 29d 50 (R n = 0.10) and t s = 19d 50 (R n = 0.01) (Fig. 12B). The thickness of the wall shear zone t s was higher by about 50% when initially loose sand was subjected to shearing (Figs. 10 and 14A) due to the smaller softening rate.
The thickness of the wall shear zone was found to be almost the same independently of pressure in the range of 50-200 kPa (Fig. 14B) due to the same rate of softening for the different pressures (Fig. 6). The thickness of shear zones increases with a decreasing rate of softening [51,56] due to the smaller post-peak stiffness. Our outcome is in contrast to DEM simulations [12,25] wherein the thickness slightly increased as the pressure decreased for R n-= 0.2-1.0 and r n = 40-100 kPa [24] and for R n = 0.02 and r n = 10-500 kPa [12]. The issue of the effect of pressure on the thickness of shear zones merits further investigations. The sphere rotations had nearly always the same sign (clockwise rotation) (Fig. 11). Only a few ones rotated in the opposite direction (Fig. 11). For R n C 0.75, the largest sphere rotation was located slightly above the bottom wall (h/d 50 = 5-6), and for R n B 0.5, it was directly located at the bottom wall (Figs. 10, 11 and 12B) where it diminished with the reduction of R n . For R n C 0.75, the sphere rotations were approaching zero at the bottom wall (the spheres were trapped in asperities). Above the wall shear zone, all spheres were almost motionless. Both the maximum horizontal displacement (Fig. 12A) and void ratio (Fig. 12C) increased in the wall shear zone with growing R n up to R n \ 0.75, respectively, and for R n C 0.75, they were almost the same. The maximum void ratio in the dilatant wall shear zone at the residual state changed between e = e 0 = 0.55 (R n = 0.01) and e = 0.85 (R n-C 0.50). The horizontal slip along the bottom obviously decreased with increasing R n (Fig. 12A), and constituted about 0% (R n = 0.75-2.0), 20% (R n = 0.50), 70% (R n-= 0.25), 93% (R n = 0.10) and 99% (R n = 0.01) of the total prescribed horizontal displacement of u x = 7 mm (Fig. 12A).
The ratio between the wall grain rotation and wall grain slip A = xd 50 /u systematically reduced with decreasing R n (Fig. 15) or increasing 1/R n (Fig. 16) The ratio xd 50 / u was approximately equal to R n at the residual state [(xd 50 )/u % R n ] and might be, thus, used as a boundary condition in micro-polar continua [55]. For the case of R n C 0.75, where all sphere rotations x tended to zero, the boundary condition could be simplified as x = 0 and u = 0 (no rotations and displacements) as in the experiment [49].  (Fig. 17) and tangential contact forces (Fig. 18) in the entire granular specimen at the residual state (front view) for the different normalized wall roughness R n and initial void ratio e o .
The distribution of both forces was similar in initially dense sand (Figs. 17A and 18A). Some differences occurred for initially loose sand only (Figs. 17B and 18B). The results of Figs. 17 and 18 evidently show that the nonuniformity of contact forces might be pronounced during wall shearing, in particular, for initially dense sand and very rough and rough walls (R n [ 0.25). The non-uniformity of contact forces increased with growing R n and diminishing e o (the distribution of contact forces for relatively smooth surfaces (R n = 0.01 and R n = 0.10) was quite uniform). The contact forces were higher at the left side wall where a passive state developed in contrast to an active state at the right side wall. In initially loose sand (R n = 1.0), the normal contact forces were more uniformly distributed along the specimen height due to higher porosity and less sphere contacts (Fig. 17B). The distribution of tangential contact forces (R n = 1.0) was similar, independently of e o (Fig. 18). Figure 19 presents a polar distribution of contact forces in the x-y plane (mean amplitude and orientation to the horizontal) at the beginning of the wall friction test after a settlement process and at the residual state. Initially, the vertical mean contact forces (with the orientation the horizontal of 90 o ) dominated due to vertical confinement r n imposed to the specimen. As a wall friction process proceeded, the direction of mean contact forces changed from a vertical to a diagonal as in other DEM results in [12,44]. Depending on the normalized surface roughness, the orientation of the mean contact forces to the horizontal at the residual state varied from 125 o (R n = 0.01) up to 160 o (R n = 0.75-2.0). For rough and very rough surfaces (Fig. 19a-e), the final maximum diagonal contact forces were 1.5 times higher than the maximum initial vertical contact forces. In initially loose sand for R n = 1.0 (Fig. 19B), due to the lower amount of interactions, the mean contact forces were greater than in initially dense sand. The final orientation angle of mean contact forces to the horizontal (e o = 0.80, R n = 1.0) was about 115 o (Fig. 19B).
The number of all sphere contacts, expressed by the coordination number c (Fig. 20), was correlated with the change of volumetric strain e v and void ratio e (Figs. 5B and 6B). For initially dense sand (R n = 1.0 and r n-= 100 kPa), the coordination number diminished with increasing of volumetric strain, corresponding to material dilatancy (Fig. 20). The maximum value of c was c = 5.25, and the residual one was c = 4.6.

Distribution of stresses and couple stresses
The stresses and couple stresses of a single sphere were calculated as [34] x c i f c j ; ð8Þ and where N-the contact points, x c i -the ith component of branch vector jointing from the centre of mass of the particle to the contact point 'c', f c j -the jth component of the total force at the contact point 'c', m c j -the jth component of the total moment at the contact point 'c' and V p -the cell volume. The values of r ij and m ij were again calculated from an averaging cell of the size 5d 50 9 5d 50 9 1d 50 being moved by d 50 . Figures 21 and 22 show the distributions of the normal stresses r ii , shear stresses r ij and stress moments m ij across the normalized specimen height h/d 50 at the specimen midpoint at residual state for u x = 7.0 mm for the different normalized wall roughness parameter R n (r 11 -horizontal normal stress, r 22 -vertical normal stress, r 12 -horizontal shear stress, r 21 -vertical shear stress, m 23 -horizontal stress moment and m 13 -vertical stress moment). The relationship between the wall shear stress-wall couple moment ratio (r 12 9 d 50 )/m 23 and the normalized wall roughness parameter R n [55] at the wall in the residual state is demonstrated in Fig. 23.
The wall shear stresses r 12 and r 21 grew with increasing wall roughness and were non-symmetric at the wall (r 12= r 12 ) (Figs. 21 and 22). The distribution of the vertical normal stress r 22 was almost the same across the specimen height in contrast to the horizontal normal stress r 11 and shear stresses r 12 and r 21 . The wall couple stresses m 13 (vertical) and m 23 (horizontal) diminished with reducing wall roughness in the range 0.25 B R n B 2.0 (for R n = 0.1 and R n = 0.01, they were low with a different sign). The distributions of stresses and couple stresses are in agreement with FEM results within micro-polar continuum [49,51,55] and other DEM results [29]. A bi-linear trend between the ratio (r 12 9 d 50 )/m 23 and R n was solely obtained for 0.25 B R n B 2.0 (Fig. 23).

Distribution of wall forces and wall moments
First, the wall forces n' and t' and wall moments m' along the grooves (Fig. 24) were calculated at the specimen midregion at the residual state for u x = 7.0 mm with the different normalized wall roughness parameter R n (e o = 0.55 and r n = 100 kPa). The wall values were computed from the area 5d 50 9 5d 50 . Next, they were transformed into a global system. Figure 25 presents the normal wall force n' in the global vertical direction, tangential wall force t' in the global horizontal direction, tangential wall moment m' in the global horizontal direction and ratio between the tangential wall force and tangential wall moment ratio (t' 9 d 50 )/m' in the global horizontal direction acting at the groove height on the granular segment. The distribution of the vertical normal wall force n' was obviously constant due to the constant vertical pressure. The relationship between the horizontal tangential wall force t' and normalized roughness parameter R n was also Normal and tangential wall forces n' and t' and tangential wall moments m' acting on inclined wall grooves bi-linear (Fig. 25Ab) as in the case of the wall friction angles u w in Fig. 8. The horizontal tangential wall force t' and horizontal tangential wall moment m' grew with increasing wall roughness up to R n = 0.75. The ratio between the horizontal tangential wall force and horizontal tangential wall moment (t' 9 d 50 )/m' was almost constant independently of R n (about 4-5) except the case of a very smooth wall with R n = 0.01 (Fig. 25c) wherein a very high value was obtained due to an extremely low value of the wall moment. The relationship (t' 9 d 50 )/m 0 = const might be used as a wall boundary condition at wall nodes in FEM using micro-polar continua.

Effect of wall groove type
Various types of a very rough wall were investigated with the normalized wall roughness parameter R n = 1.0 (Fig. 26). The effects of the grooves' distance s g (varying between 2 9 d 50 and 16 9 d 50 ) and grooves' inclination to the bottom a g (45 o or 27 o ) (Fig. 3) were analysed in detail (Figs. 26, 27 and 28).
The peak wall friction angle, maximum volumetric strain and thickness of the wall shear zone decreased with increasing s g and with decreasing a g (Figs. 27 and 28). The effect of the groove distance s g proved to be small and the effect of the groove inclination a g was strong on the peak wall friction angle (Fig. 27). The effect of the groove Fig. 25 Relationships between normal wall force n' in vertical direction (a), tangential wall force t' in horizontal direction (b) (A), tangential wall moment m' in horizontal direction (B) and ratio between tangential wall force and tangential wall moment (t' 9 d 50 )/m' in horizontal direction versus normalized roughness parameter R n at grooves' height (C) at residual state (e o = 0.55 and r n = 100 kPa) distance started to be visible for s g C 89d 50 . The peak wall friction angle was lower by 5% for s g C 89d 50 than for s g = (2-4) 9 d 50 . The asperity angle of a g = 45 o yielded the higher maximum wall friction angle than a g = 27 o by 10%. The residual wall friction angle was not affected by the grooves' type. Since the effect of the groove inclination was stronger than the groove spacing, the modelling of the wall roughness with the same groove inclination but the different groove distance is more reasonable [25,63]. In the next calculation step, the regular triangular grooves with a different asperity height will possess the same inclination angle as in [25]. The DEM wall friction simulations will be also performed with sand grains in the form of non-symmetric clumps [28].

Conclusions
The results of a series of DEM simulations of varying interface roughness were presented for cohesionless sand. Comparisons of results between the DEM simulations and previously reported physical tests for different interface roughness showed good qualitative agreement in measured wall friction angles. Some main conclusions can be offered from our simulations: -The DEM model produced trends consistent with those observed in physical experiments. The normalized interface roughness had a huge influence on the mobilized wall friction angle and thickness of the wall shear zone.
-The wall friction resistance increased with increasing wall roughness. The peak and the residual wall friction angle rapidly increased with increasing roughness parameter up to a particular value of the normalized wall roughness (R n = 0.75). The relationship between the peak/ residual wall friction angle and normalized wall roughness was bi-linear. The peak wall friction angle increased with decreasing normal stress and increasing initial void ratio. The residual wall friction angle was the same independently of the initial void ratio and pressure.
-The wall friction resistance was strictly combined with the thickness of the wall shear zone wherein pronounced grain rotations occurred. The wall friction resistance grew with the increasing normalized wall roughness parameter up to R n = 0.75. The wall shear zone extent expanded in a bi-linear way from 1 9 d 50 up to 14 9 d 50 for R n = 0.01-2.0. The wall shear zone thickness increased with growing initial void ratio and was almost the same for the different vertical pressure (R n = 1).
-The effect of the groove distance s g proved to be small and the effect of the groove inclination a g proved to be significant on the peak wall friction angle and the thickness of the wall shear zone. The maximum interface efficiency was achieved in DEM studies for the asperity spacing to the mean grain diameter ratio between 1.0 and 4.0, and the asperity height to mean grain diameter ratio equal to or greater than 0.75. The asperity angle of 45 o yielded the highest efficiency.
-For very rough walls (R n C 0.75), the largest sphere rotation was located slightly above the bottom wall (h/ d 50 = 5-6), and for R n B 0.5, it was directly located at the bottom wall where it diminished with the reduction of R n . For R n C 0.75, the sphere rotations were approaching zero at the wall (the spheres were trapped in asperities).
-At the residual state, the ratio between the wall grain rotation and wall grain slip xd 50 /u was found to be directly linked with the normalized wall roughness parameter R n . The ratio between the tangential wall force and wall moment (t' 9 d 50 )/m' along the horizontal wall was almost constant independently of R n .
-The non-uniformity of contact forces in granular specimens increased with the growing wall roughness and decreasing initial void ratio. The orientation of mean contact forces at the residual state also grew with the increasing wall roughness (from 125 o up to 160 o for Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.