Comparative 3D DEM simulations of sand–structure interfaces with similarly shaped clumps versus spheres with contact moments

Three-dimensional simulations of a monotonic quasi-static interface behaviour between initially dense cohesionless sand and a rigid wall of different roughness during tests in a parallelly guided direct shear test under constant normal stress are presented. Numerical modelling was carried out by the discrete element method (DEM) using clumps in the form of convex non-symmetric irregularly shaped grains. The clumps had an aspect ratio of 1.5. A regular grid of triangular grooves (asperities) along the wall with a different height at the same distance was assumed. The numerical results with clumps were directly compared under the same conditions with our earlier DEM simulations using pure spheres with contact moments with respect to the peak and residual interface friction angle, width of the interface shear zone, ratio between grain slips and grain rotations, distribution of contact forces and stresses. The difference between the behaviour of clumps and pure spheres with contact moments proved to be noticeable in the post-peak regime due to a different particle shape. The rolling resistance model with pure spheres was proved to be limited for capturing particle shape effects. Three different boundary conditions along the interface were proposed for micropolar continua, considering grain rotations and grain slips, wall grain moments and wall grain forces, and normalized interface roughness. The numerical results in this paper offer a better understanding of the interface behaviour of granular bodies in DEM and FEM simulations.


Introduction
The safety of engineering systems adjacent to soils in static and dynamic conditions is strongly affected by the strength of soil-structure interfaces. The behaviour of interfaces is characterized by the formation of shear zones with a certain thickness in the soil adjacent to structures, i.e. thin zones of intense shearing with both pronounced grain rotations and volume changes [52,58]. The thickness of wall shear zones determines the forces transferred from soils to structures. To properly describe the width of shear zones in finite element simulations, continuum models must be enriched by a characteristic length of micro-structure [56][57][58][59] by means of, e.g., micropolar, second-gradient, non-local or viscous terms (e.g. [2,40,54,68]). Micropolar models are, in particular, physically realistic for cohesionless soils since they include grain rotations that are visible in experiments [52]. In those models, the characteristic length of micro-structure directly corresponds to a mean grain radius/diameter (e.g. [56][57][58][59]). In the case of enriched continua including, e.g. second-gradient, non-local or viscous terms, the characteristic length is determined with an inverse identification process of experimental data [2], [54]. Up to now, there do not exist consistent micropolar continuum boundary conditions for walls of different roughness [10,11,20,58,60]. There are some elastoplastic and hypoplastic models in the classical continuum for describing granular soil-structure interface behaviour (e.g. [37,45,49,51]). However, Coulomb's friction law is usually used in problems of the soil-structure interface, based on the assumption of a constant ratio between the shear and normal stresses on the interface. This assumption is incorrect since: (1) the mobilized friction angle can change significantly during shearing along rough or very rough walls in contact [52,58] and (2) the mobilized wall friction angle is not a state constant since it depends on boundary conditions, contact pressure level, initial stress state and specimen size [52,58]. Therefore, the laboratory wall friction angles obtained from different tests (in particular, their peak values) cannot always be directly transferred to other boundary value problems [52].
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 [9]. An asperity angle of 50°or greater yielded the maximum efficiency for any given asperity spacing or height [9]. The maximum interface friction angle was found to be larger than the internal friction angle of soils due to the passive resistance caused by surface asperities [21,22,52]. It varied in a bi-linear fashion as a function of the normalized roughness [19,42,50]. The critical normalized roughness was mainly about 0.40-0.50 and above this value, the maximum interface friction angle insignificantly increased. The literature overview on the interface behaviour between cohesionless sand and rigid wall of different roughness with respect to numerical solutions within continuum mechanics and discrete mechanics was given in [15]. There exist many DEM studies of sand-structure interfaces under 2D (e.g. [13,14,16,22,64,65,74]) and 3D conditions [12,15,24,69,73]. In those studies, disks with contact moments were used in 2D simulations [13,14,16,22,64,65,74] or spheres with contact moments were used in 3D simulations [12,15,24,69]. The clumps were either assumed in 2D [22] or 3D simulations [73]. Zhou et al. [73] performed comprehensive 3D simulations of the soil-structure interface behaviour with different interface roughness, particle sphericity and initial fabric, using three types of clustered particles. The shear strength of the interface was found to grow as the interface roughness and particle orientation increased and particle sphericity reduced. However, the DEM studies directly comparing the behaviour of spheres including contact moments with clumps in the soil-structure interface problems under the same conditions are missing. Those comparative studies with a different particle shape are needed since the majority of DEM simulations are carried out using disks/spheres with contact moments to strongly shorten the computation time. This approach is physically less realistic than the approach with clumps. A comprehensive comparison between the rolling resistance and particle shape effects was carried out with DEM for triaxial compression tests only, without taking shear localization into account [31,71]. It was shown that the rolling resistance model was limited in its ability to capture particle shape effects [71]. The existing DEM simulations also exhibited a bi-linear relationship between the interface resistance and normalized interface roughness [15,24,69,73]. The critical normalized roughness was found to be in 3D studies about 0.40-0.75 [15,24,69,73].
In our earlier paper [15], the monotonic quasi-static interface behaviour between cohesionless sand and a rigid wall of different surface topography using a direct shear box under conditions of constant normal stress was comprehensively analysed using 3D DEM (the particle-based open-source code YADE). To approximately simulate the irregularity of sand particles' shape of 'Karlsruhe' sand, spheres with contact moments were used. The artificial wall surface topography was created by regularly arranged triangular grooves (asperities) in the form of a standard saw-tooth surface. Some calculation results were compared with our corresponding monotonic quasi-static sandstructure interface tests on cohesionless 'Karlsruhe' sand in a parallelly guided direct shear box, performed at Karlsruhe University [52,58]. They were also compared with reported experimental results for similar interface roughness in the literature. The calculations were carried out with the different interface roughness, expressed by the varying height, spacing and inclination of wall grooves, initial void ratio of sand and vertical pressure. The effect of the groove distance was found to be small up to a certain value and the effect of the groove height and groove inclination was significant on the peak interface friction angle and the thickness of the shear zone. In addition, three wall boundary conditions for 2D micropolar continua were proposed: (1) the ratio between the micropolar rotation and slip, (2) the ratio between the interface shear stress and interface couple moment and (3) the ratio between the tangential interface force and tangential wall moment.
In the current paper, the 3D interface behaviour was investigated under the same conditions, using both the same direct shear box and granular material. However, instead of pure spheres with contact moments [15], we assumed irregular grains of the same shape composed of sphere clumps that more realistically reproduced the real grain shape. The aspect ratio of grains was 1.5. In the present calculations, the initial void ratio and vertical pressure were also the same. The artificial wall surface topography was again created by regularly arranged triangular grooves (asperities) at the same spacing [15]. The attention was focused on differences in the sand-structure interface behaviour, described by two different DEM approaches, with respect to the evolution and magnitude of the mobilized wall friction angle, wall shear zone, interface displacements, interface rotations, stresses, interface forces/moments and contact forces in the granular assembly for the interface of the different roughness. The goal was to find if a rolling resistance model with spheres is able to capture particle shape effects. In addition, interface boundary conditions of different roughness were assumed for clumps within micropolar continuum by considering grain rotations, slips, moments and forces. They were directly connected with the normalized interface roughness, using simple relationships. As compared to [15], the ratio between the slip and the total imposed displacement along the interface was proposed instead of the ratio between the horizontal shear stress multiplied by the mean grain diameter and the horizontal couple stress. The ratios between interface grain rotations and interface grain slips and between interface grain forces and interface grain moments were not measured in experiments yet. The particle slip was separately measured in experiments by DeJong et al. [6].
The paper includes the following novel points: (1) comprehensive 3D analyses of the interface behaviour with the real mean grain diameter of the sand using irregularly shaped grains by taking into account the effect of different interface roughness, (2) a detailed comparison of the interface behaviour between grain clumps of the same shape and pure spheres with contact moments under the same conditions to show the particle shape effect on the shear resistance and shear localization (such a comparison was not performed for sand-structure interface problems yet) and (3) a proposal of simplified interface boundary conditions of different roughness for micropolar continua by considering grain rotations, slips, moments and forces related to the normalized interface roughness. The findings presented in this paper can help to better understand both: (1) the sand-structure interface mechanism when using two different DEM approaches and (2) the boundary conditions along the interface of different roughness for micropolar continua to simulate large soil-structure systems with the finite element method.
The paper is arranged as follows. In Sect. 2, the proposed numerical DEM framework was presented in detail. The model calibration is shown in Sect. 3. Results of the effect of the interface roughness on friction angle and shear localization with key findings are described in Sect. 4. Section 5 discusses interface boundary conditions in the frame of the micropolar (Cosserat) continuum. The main results were summarized and some conclusions are stated in Sect. 6.

3D DEM model
DEM is a very efficient and powerful numerical tool for modelling granular assemblies. The 3D spherical discrete element model YADE developed at the University of Grenoble [27,28,47] was employed to simulate the interface behaviour. To simulate the irregularity of the particles' shape of 'Karlsruhe' sand, clumps composed of spheres were used [30]. YADE takes advantage of the socalled 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 mechanical response of the contact model when using clumps is presented in Fig. 1. The DEM model for spheres with contact moments was presented in [15]. Below, the basic equations of the DEM model for clustered particles are summarized [27,28,47]: 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, The three 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). In addition, a particle radius R, particle mass density q and numerical damping parameter a d are required. The DEM material parameters: E c , m c and l were calibrated using the corresponding homogeneous axisymmetric triaxial laboratory test results on 'Karlsruhe' sand with the different initial void ratio and lateral pressure [31] that were compared with the corresponding experiments [26,66]. The procedure for determining the material parameters in DEM was described in detail by Kozicki et al. [30,31]. Note that the representative elastic contact moduli E c and m c are different from the elastic moduli of grains. The effect of damping (if a d-B 0.10) is negligible in quasi-static calculations [30,31].
All DEM calculations on the interface behaviour were performed with convex non-symmetric irregularly shaped grains of a diameter d clump , composed of 4 spheres of different diameters (d sphere = 0.50-0.60 9 d clump ) (Fig. 2). Each granular assembly was prepared by putting clumps of a random diameter according to the grain distribution curve. The aspect index (the ratio between the maximum and minimum clump diameter) was 1.50, the convexity index '1' (the ratio between the smallest sphere volume encompassing the cluster and the cluster volume) was 2.07 and the convexity index '2' (the ratio between the smallest convex volume encompassing the cluster and the cluster volume) was 1.16.

Model calibration based on triaxial compression
To find the material parameters in DEM (E c , m c and l), a series of numerical homogeneous quasi-static triaxial compression tests with rigid smooth walls on cohesionless sand were initially performed [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,66] (Fig. 3). The triaxial compression test is the most homogeneous laboratory test in soil mechanics and hence is frequently used for soil calibration. The index properties of 'Karlsruhe' sand are: the 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 c d max = 17.4 kN/m 3 , minimum void ratio e min = 0.53, minimum specific weight c d min = 14.6 kN/m 3 and maximum void ratio e max = 0.84. The sand grains were classified as sub-rounded/sub-angular. The material parameters could not be calibrated with our interface shear tests [58] since we did not have full knowledge of the wall surface whose roughness was very non-uniform and random.
For DEM simulations, a cubical specimen of size 10 9 10 9 10 cm 3 (composed of about 8000 grain clumps as in [15]) (Fig. 3a). The grain diameter of sand linearly varied between 2.5 and 7.5 mm and its mean grain diameter was d 50 = 5 mm (10 times larger than the real one d 50 = 0.5 mm). The mass density was 2600 kg/m 3 . The DEM simulations were carried out for initially very dense sand (initial void ratio e o = e min = 0.53) with three different lateral pressures r c = 50 kPa, 200 kPa and 500 kPa (Fig. 3b, c) assuming the material parameters listed in Table 1. A satisfactory agreement was obtained, in particular, for the stress-strain curves (Fig. 3b). This agreement with laboratory tests was slightly worse than for pure spheres with contact moments [15] since the DEM model for spheres had a higher number of material parameters [15] (five against three). When using clustered particles, a better agreement with experiments is solely possible with more realistic grain shapes, obtained with the aid of, e.g. 3D micro-CT images [25,41]. The effect of E c and l on the sand behaviour is presented for r c = 200 MPa in Fig. 4. The lower the parameter E c , the lower are the global elastic modulus E, maximum internal fiction angle / max and volumetric strain e v and the higher is the vertical normal strain e 1 corresponding to / max . The reduction of E c by factor 3 diminishes E by 60%, / max by 3% and e v by 10%, and increases e 1 corresponding to / max by 50%. The lower the parameter l, the lower are the maximum internal fiction angle / max and volumetric strain e v . The reduction of l by 20% diminishes / max by 10% and e v by 20%. The residual internal fiction angle / res is not affected by E c and l. The vertical normal strain e 1 corresponding to / max is 5-7% (Fig. 4A).
A simple linear elastic normal contact was chosen to capture on average various contact possibilities present in real sands [39] and to reduce the computation time [67]. We performed comparative simulations for spheres with  contact moments between a linear and nonlinear contact model [36] following Hertz [18] and Mindlin and Deresiewicz [35], based on a plane compression test without shear localization. The comparative analyses were carried out only with respect to macroscopic quantities (measured in real laboratory tests) due to the lack of experimental data on microscopic quantities (Fig. 5). The results of triaxial compression indicated some negligible differences between the two contact models for spheres and contact moments. Both the global maximum mobilized internal friction angle and dilatancy angle for initially dense sand were almost the same using both approaches, and the global macroscopic elastic parameters (elastic modulus and Poisson's ratio) differed about by 10% only. For initially dense sand and confining pressure r c = 200 kPa, the global macroscopic elastic parameters were E = 80 MPa and m = 0.25 using spheres with contact moments and a linear contact model, and E = 70 MPa and m = 0.22 using spheres with contact moments and a nonlinear contact model. The global maximum mobilized internal friction angle u max and dilatancy angle were u max = 42°and w = 30°(spheres with contact moments and a linear contact model), and u max = 42.1°a nd w = 30°(spheres with contact moments and a nonlinear contact model) (Fig. 5). Although a nonlinear contact model is more realistic (in particular for small loads) [39], a linear elastic contact law provides similar DEM results with the significantly reduced computation time [67] and therefore was used in current simulations. The elastic constants of the grain contact in our DEM model do not correspond to the elastic constants of the spheres' material in the Hertz law (the parameter K n in Eq. 3 is several times larger than the mean normal stiffness of the spherical grain material) [1].

DEM model
The 3D granular specimen included 80 000 grain clumps (Fig. 6). The specimen length (l = 100 mm) and the height (h = 20 mm) (Fig. 6A) were the same as in the experiment [52,58]. The width of the granular specimen with clumps was D = 5 mm (10 9 d 50 ) in all DEM simulations instead of D = 100 mm as in the experiment, based on preliminary simulations of a direct interface shearing test using pure spheres with contact moments [15]. The reduction of the specimen width from D = 100 mm down to D = 5 mm did not affect the results [15] but strongly diminished the computation time. The gap equal to the maximum grain diameter was left between the bottom and vertical walls (see the zoom in Fig. 6A) as in the experiment [52,58] to prevent locking of particles at the bottom corners during shearing. The sand leakage during shearing had a minor effect on void ratio and volumetric strain (\ 0.5%). The maximum number of clumps beyond the box after the test was solely 350. The comparative calculations without a gap indicated similar results; however, the evolution of the interface friction angle showed more fluctuations caused by grain interlocking at ends. The grain diameter of sand linearly varied between 0.25 mm (d min = 0.25 mm) and 0.75 mm (d max = 0.75 mm). The mean grain diameter of d 50 = 0.5 mm was as in the experiment [52]. The sand was horizontally sheared along the bottom wall under a constant velocity of the shear box du x /dt = 0.5 mm/s (from the left to the right, u x -the positive horizontal displacement) (Fig. 6). 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 ) [34]. The rigid bottom wall was fixed. The top and side rigid walls were movable and smooth. The frictional coefficient along those walls was tanu w = tanl = tan18°= 0.32. In DEM simulations, the interface roughness was simulated in a various way [12,23,64]. As in [15], a regular saw-tooth surface with varied asperity height and the same distance was chosen. The regularly arranged triangular grooves (asperities) had the same distance but a different inclination. The interface roughness was characterized by the normalized interface roughness parameter R n = h g /d 50 , where h g is the groove height and d 50 denotes the mean grain diameter [58,61,69]. The parameter R n was 2.0, 1.0, 0.75, 0.50, 0.25, 0.10 and 0.01 (Fig. 6B). The groove distance s g was always the same (s g = 2 9 d 50 ). The groove inclination to the bottom a g diminished with decreasing R n (e.g. a g = 45°for R n = 1.0) (Fig. 6B). The constant uniform vertical pressure r n was applied to the top area of the shear box. The interface friction angle between particles and grooves/walls was assumed to be u w = 18°. The comparative calculations with particle clumps versus pure spheres with contact moments [15] were limited in the current paper to initially dense sand with e 0 = 0.55 and one vertical pressure of r n = 100 kPa (as in [15]). The effect of the different vertical pressure r n and initial void ratio e 0 on the sand behaviour was investigated for spheres with contact moments in [15]. For clumps, a similar effect is expected. The same initial void ratio was assumed in both approaches since the calculations concerned the same sand. We did not have problems generating numerically this initial void ratio for clumps due to an overlap of particles. The relative initial density of the granular material should be controlled if different sands are compared with each other [38]. The relative density was controlled in DEM simulations of the interface shear behaviour for different sands by, e.g. Zhou et al. [73]. The interface friction angle u was calculated as u = arctan(T/N) (T-the sum of horizontal forces acting on the top and lateral walls and N-the sum of vertical forces acting on the top wall). The calculation time for clumps was larger by factor 4 as compared to that for spheres with contact moments.

Effect of interface roughness on mobilized friction angle
The normalized interface roughness parameter R n varied between R n = 0.01 and R n = 2.0 for e 0 = 0.55 and r n-= 100 kPa (Fig. 7) as in DEM analyses with pure spheres including contact moments [15]. The values of the peak interface friction angle u w,max , residual interface friction angle u w,res and volumetric strain e v are given in Table 2 as compared to pure spheres with contact moments [15]. The relationship between the computed values of u w,max and u w,res and R n and grain shape is shown in Fig. 8.

Evolution of mobilized interface friction angle and volumetric strain
The evolution shape of the curves u w = f(u x ) and e v = f(u x ) ( Fig. 7) for the initially dense sand are similar as for pure spheres with contact moments [15]. The peak u w,max and the residual interface friction angle u w,res continuously grew with increasing roughness parameter. The peak values of u w,max were very similar for both clumps and pure spheres (Table 2); however, the residual values of u w,res and residual volumetric strains were slightly lower for clumps with the normalized interface roughness parameter R n C 0.25 (e.g. they were lower by 3°for R n C 0.5 and by 6°for R n = 0.25). The clear critical surface roughness R n(crit) concerning the peak interface friction angle u w,max was not reached with clumps in contrast to pure spheres with contact moments wherein the clear critical surface roughness was obtained for R n(crit) = 0.50-0.75, beyond which their effect became negligible ( Table 2). The relationship between the values of u w,max and R n was more parabolic for clumps in opposite to a bi-linear one for pure spheres although the changes of the peak interface friction angle were small for clumps above the value of R n C 0.75 (Fig. 8). The relationship between the values of u w,res and R n was, however,  Relationship between peak interface friction angle u w,max and residual interface friction angle u w,res and normalized interface roughness parameter R n from DEM for: a) clumps and b) pure spheres with contact moments [15] (black colour indicates u w,max and red colour indicates u w,res ) (color figure online) 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 interface roughness parameter R n using clumps of Fig. 2 (values for pure spheres with contact moments [15]  bi-linear for both clumps and pure spheres [15]. The value of u w,max = 188 for R n = 0.01 was obviously equal to the assumed interface friction angle between particles and grooves. The slightly lower values of R n(crit) were obtained in the experiments by Hu and Pu [19,24] (R n(crit) = 0.40) and DEM calculations by Jing et al. [23,73], Zhou et al. [69,73] (R n(crit) = 0.375-0.50) and Zhang and Evans [69] (R n(crit) = 0.40). A similar value of R n(crit) was obtained in the laboratory tests by Su et al. [50] (it was between 0.5 and 1.0). The volumetric strain of sand also changed with increasing R n (Fig. 7B); its greatest value at the residual state was lower by 20% than for pure spheres with contact moments [15]. Figure 9 shows the effect of different frictional angles between sand particles and lateral walls of the shear box (10°and 26°) on the DEM results. The effect was insignificant on the residual friction angle at the interface and volume changes of the granulate. The difference between the maximum interface friction angles was 5% only.
As in [15], the DEM results were also compared with our interface experiments [52,58]. In experiments, the interface roughness was classified as smooth (0 \ R n-B 0.1 9 d 50 ), rough (0.1 9 d 50 \ R n \ 0.5 9 d 50 ) and very rough (R n C d 50 ). The rough interface was obtained in a corrosion chamber, and a very rough interface 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. Thus, the experimental interface roughness was different. The DEM calculation results were compared with the experiments in Fig. 10. A satisfactory agreement with the experimental results was obtained with respect to the shape of curves u w = f(u x ) and e v = f(u x ) and peak interface friction angles u w,max . However, the calculated residual interface friction angles u w,res were much smaller than the experimental values for the rough and very rough interface. The calculated displacements corresponding to u w,max were also smaller (by factor 3) for the rough and very rough interface than in the experiment. Thus, the numerical initial sand response was too stiff. The calculated volumetric strain e v was also too large by factor 1.5 as compared to the experiments for the rough and very rough wall. Those discrepancies were caused by the fact that the experimental interface roughness and particle shapes were not faithfully reproduced in DEM simulations. A drop/growth of the parameter l causes a decrease/increase of the parameters u w,max and e v , whereas the value of u w,res remains the same (Fig. 4A). A decrease of E c improves the initial material response, volumetric strain and the location of u w,max but worsens u w,max (Fig. 4) as compared to the experiments (Fig. 10) without affecting again the value of u w,res . Hence, the parameter l must be higher. On the other hand, a small value of E c is not realistic with respect to triaxial compression laboratory tests (Fig. 4). To increase u w,res for better fitting with the experiments (Fig. 10), the grain shape with a higher aspect ratio should be used [30]. Note that the experimental results in the initial shear phase may be affected by some test inaccuracies that contribute to a weaker sand response at the beginning of deformation. The measured horizontal displacement u x corresponding to u w,max for very rough walls was about 1.5 mm (Fig. 10). The sand strain corresponding to this displacement, approximately calculated as e x = u x /t s= 1.5 mm/6 mm = 0.25 (t s the thickness of the shear zone at the interface, Fig. 13), is much too high by comparing it with triaxial laboratory test results (Fig. 3).

Distribution of particle rotation, particle displacement and void ratio
The effect of the normalized interface 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. 11. The values of rotations x were calculated from a cubic averaging cell of the size 5d 50 9 5d 50 9 5d 50 being moved by d 50 . In the averaging cell, the centres of spheres were considered independently of their diameter (clockwise rotation is positive). Figure 12 presents the zoom on single sphere rotations in the granular segment at the front side (mid-length). In Fig. 13, the relationship between the shear zone thickness t s and normalized interface roughness parameter R n is demonstrated as compared to pure spheres with contact moments [15]. Figure 14 shows the distribution of sphere horizontal displacements u x , 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 . The relationship between the horizontal interface slip and total prescribed horizontal displacement and normalized interface roughness parameter R n is presented in Fig. 15. Figure 16 demonstrates the evolution of the ratio between the interface grain rotation and interface grain slip A = (xd 50 )/u during interface shearing at the box mid-length. The distribution of the residual ratio between the grain rotation multiplied by the mean grain diameter xd 50 and the sphere slip u (xd 50 )/u across the specimen height at the specimen midpoint for the different parameter R n is shown in Fig. 17.
The ratio A at the residual state is shown in Table 3.
The particle rotations and increasing void ratio are the best indicators for shear localization [15,52,[56][57][58][59]. Based on particle rotations (Fig. 11), it can be seen that during interface shearing with the different normalized interface roughness parameter R n , an almost horizontal dilatant shear zone was created at the interface (R n [ 0.01). The grain rotations had nearly always the same positive sign (clockwise rotation) (Fig. 11). More grains rotated in the opposite direction (Fig. 12) as compared to pure spheres with contact moments [15]. The thickness of the shear zone was estimated, based on an inflection point in the distribution of sphere rotations x where the rotation was x B 5% of the maximum particle rotation x max in the wall shear zone (Fig. 14B). The thickness of the shear zone t s at the interface increased in an approximate parabolic way (Fig. 13a). It was equal to t s = 12 9 d 50 (R n = 2.0), t s = 11.0 9 d 50 (R n-= 1.0), t s = 10 9 d 50 (R n = 0.75), t s = 8 9 d 50 (R n = 0.50), t s = 6 9 d 50 (R n = 0.25), t s = 3 9 d 50 (R n = 0.10) and t s-= 1 9 d 50 (R n = 0.01) (Fig. 11) (for pure spheres, it changed in a bi-linear one, Fig. 13b). The calculated maximum thickness was in agreement with the maximum experimental value of 9 9 d 50 [6]. The thickness was lower than for spheres [15] by 15-40% for R n C 0.50 (Fig. 13) due to higher grain rotations at the wall (Fig. 14B). Similarly, as for pure spheres with contact moments [15], the largest grain rotation was located slightly above the interface (h/ d 50 = 3-5) for R n [ 0.5, and for R n B 0.5 it was directly located at the interface (Fig. 14B) where it diminished with the reduction of R n . For R n [ 0.5, the sphere rotations approached zero at the bottom wall (equal to the situation with particles trapped in asperities). Above the shear zone, all spheres were almost motionless (Figs. 11 and 14B). The maximum residual grain rotation for R n C 0.75 was equal to about 60°, whereas for pure spheres with contact moments, it was about 45° [15].
Both the maximum horizontal displacement (Fig. 14A) and void ratio (Fig. 14C) increased in the shear zone at the interface with growing R n . The distribution of the horizontal displacement and void ratio across the specimen height was similar as in interface shear laboratory tests using the particle image velocimetry PIV [5][6][7]. The maximum void ratio in the dilatant wall shear zone at the residual state changed between e = 0.55 (R n = 0.01) and e = 0.86 (R n = 2.0). The horizontal slip u along the wall obviously decreased with increasing R n , and was about 0% (R n = 2.0), 7% (R n = 1.0), 14% (R n = 0.75), 35% (R n-= 0.50), 70% (R n = 0.25), 95% (R n = 0.10) and 99% (R n-= 0.01) of the total prescribed horizontal displacement of u x = 7 mm (Fig. 15a). The ratio B = u/u x was similar (R n-B 0.25) or slightly higher (R n [ 0.25) than for pure spheres with contact moments (Fig. 15b). The evolution of ratio B was hyperbolic in the range 0 B R n B 2.0 and almost linear in the range 0 B R n B 1.0 as in interface shear laboratory tests using the particle image velocimetry PIV [6].
The ratio between the interface grain rotation multiplied by the mean grain particle and interface grain slip A = xd 50 / u was almost the same during shearing (Fig. 16). The ratio A = xd 50 /u at the residual state systematically reduced with decreasing R n ( Table 3 50 )/u % R n ] for R n B 0.75 [15]. It was caused by greater grain rotations of clumps than of pure spheres with contact moments. Thus, the ratio A also depends on the grain shape. Figure 18 demonstrates the 3D distribution of the normal contact forces in the entire granular specimen at the residual state (front view) for the different normalized interface roughness R n .

Distribution of contact forces
The normal contact forces were higher at the left side wall where a passive state developed in contrast to an active state at the right lateral wall. The non-uniformity of contact forces increased with growing normalized interface roughness R n . It was pronounced during shearing at the very rough and rough interface (R n C 0.50). It was higher for pure spheres with contact moments (Fig. 18) [15]. Figure 19 presents a polar distribution of contact forces in the x-y plane (mean amplitude and orientation to the horizontal) [24,46] at the beginning of the test after a contracting process (settlement) and at the residual state for clumps. Initially, the vertical mean contact forces (with the orientation to the horizontal of 90°) dominated due to vertical pressure r n imposed on the specimen. As a sandstructure interface process proceeded, the direction of mean contact forces changed from a vertical to a diagonal as in other DEM results in [15,24,46]. Depending on the normalized interface roughness, the orientation of the mean contact forces to the horizontal at the residual state varied from 135°(R n = 0.01) up to 140°(R n = 0.25-2.0) for clumps and from 125°(R n = 0.01) up to 160°(R n-= 0.75-2.0) for pure spheres with contact moments [15]. For rough and very rough surfaces, the final maximum Fig. 13 Relationship between shear zone thickness t s and normalized interface roughness parameter R n at residual state from DEM (e o = 0.55 and r n = 100 kPa) for: a) clumps of Fig. 2 and b) pure spheres with contact moments [15]   diagonal contact forces were twice higher than the maximum initial vertical contact forces. The similar visible differences between clumps and spheres with contact moments regarding induced anisotropy of contact normal forces were also indicated in [71] for a triaxial compression test. They grew with an increase in the aspect ratio of clumps.
The number of all grain contacts, expressed by the coordination number c in initially dense sand (R n = 1.0 and r n = 100 kPa), is shown in Fig. 20 for clumps (Fig. 20a) and pure spheres with contact moments (Fig. 20b) [15]. The coordination number is an important measure to quantify the internal structure of granular materials.
The coordination number c was higher (by 55-70%) in clumps than in pure spheres with contact moments (Fig. 20). The maximum value of c was c = 9.0 (clumps) and c = 5.25 (pure spheres) and the residual one was c = 7.2 (clumps) and c = 4.6 (pure spheres). The coordination number diminished with shear due to material dilatancy, reaching an asymptote at the residual state (Fig. 7B).

Distribution of stresses and couple stresses
The stresses of a single sphere were calculated as [33] where N is the number of 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' and V p the cell volume. The values of r 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 distribution of the normal stresses r ii and shear stresses r ij across the normalized specimen height h/d 50 at the specimen mid-point at residual state for the horizontal displacement u x = 7.0 mm for the different normalized interface roughness parameter R n (r 11 horizontal normal stress, r 22 vertical normal stress, r 12 horizontal shear stress, r 21 vertical shear stress).
The stress distribution shapes were similar to those for pure spheres with contact moments [15]. The interface shear stresses r 12 and r 21 grew with increasing interface roughness and were slightly non-symmetric at the interface (r 12 = r 12 ) (Fig. 21). The non-symmetry was higher by 5-10% than for pure spheres [15]. 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 calculated stress distributions are in agreement with FEM results within micropolar continuum [52,55,60] and with other DEM results [29].

Distribution of grain forces at interface
First, the interface grain forces n and t and interface grain moments m along the grooves (Fig. 23) were calculated at the specimen mid-region at the residual state for u x-= 7.0 mm with the different normalized interface roughness parameter R n . The moments were evaluated by multiplying the tangential and normal contact forces by their distance from the gravity centres r t and r n of clumps (Fig. 23). The values were computed from the area 5d 50 9 5d 50 . Next, the forces and moments were transformed into a global system. Figure 24 presents the normal interface force n 0 in the global vertical direction, tangential interface force t 0 in the global horizontal direction, tangential interface moment m 0 in the global horizontal direction and ratio between the tangential interface force    Fig. 2 and b) pure spheres with contact moments [15] (red colour corresponds to normal contact forces higher than mean value, maximum value of forces is 0.50 N) (color figure online) and tangential interface moment ratio C = (t 0 9 d 50 )/m 0 in the global horizontal direction at the groove height in the granular segment. The distribution of the vertical normal wall force n 0 was obviously constant due to the constant vertical pressure (Fig. 24Aa). The relationship between the horizontal tangential interface force t 0 and normalized roughness parameter R n was also bi-linear (Fig. 24Ab) as in the case of the residual interface friction angles u w.res in Fig. 8. The horizontal tangential interface moment m 0 grew almost parabolically with increasing interface roughness 0 B R n-B 2.0 or linearly with increasing interface roughness 0 B R n B 1.0 (Fig. 24B). The evolution of the ratio between the horizontal tangential interface force multiplied by the mean grain diameter force and horizontal tangential  interface moment, C = (t 0 9 d 50 )/m 0 , was similar to this for pure spheres with contact moments [15]. It changed in a hyperbolic way between 1.0 (R n = 2) and 40 (R n = 0.01) and between 1.0 (R n = 2) and 8.0 (R n = 0.10) (Fig. 24C). In the case of pure spheres with contact moments, it varied also in a hyperbolic way between 4.5 (R n = 2) and 40 (R n = 0.01) and between 4.5 (R n = 2) and 7.0 (R n = 0.10) [15].

Micropolar boundary conditions
Modelling of shear zones within continuum mechanics using FEM can be only performed with constitutive models enriched by a characteristic length of micro-structure [52,56,57,59]. 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 [55]. 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 [20,52,58], strain gradient elasto-plasticity [63] and micropolar hypoplasticity [53]. Some micropolar boundary conditions were proposed in [10,11,20,60] for describing the interface roughness. In [60], the 2D boundary conditions along the horizontal rigid interface included two ratios connected to the normalized interface roughness (a ratio of the micropolar rotation multiplied by the mean grain diameter and the horizontal displacement and a ratio between the horizontal shear stress multiplied by the mean grain diameter and the horizontal couple stress). In [10,11,20], as a boundary condition in 2D analyses, the ratio between the slip and the total imposed displacement at the interface was assumed instead of the ratio between the horizontal shear stress multiplied by the mean grain diameter and the horizontal couple stress.
Based on present DEM simulation results, the following three boundary conditions at the interface, based on the ratios A = (xd 50 )/u (Fig. 17, Table 3), B = u/u x (Fig. 15) and C = (t 0 9 d 50 )/m 0 (Fig. 24C) may be used in 2D conditions. Those ratios slightly depend on the particle shape. For the case of R n C 0.75, where all grain rotations x tend to zero at the interface, the boundary condition can be simplified as x = 0, u = 0 and v = 0 (no rotations and displacements) as in the experiment [52,58].
The calculated ratios A (Fig. 17, Table 3), B (Fig. 15) and C (Fig. 24C) for clumps of Fig. 2 may be approximated yet in the range 0 B R n \ 0.75 as All ratios depend on a simple way on R n only. The proposed boundary conditions may be prescribed in nodes of FE meshes along the interface in micropolar simulations (the particle rotation x may be replaced by the micropolar (Cosserat) rotation x c [29,58,60]). The boundary conditions will be also checked in DEM interface simulations by assuming different grain shapes, grain size distributions and mean grain sizes. Relationships for clumps between normal interface force n 0 in vertical direction (a), tangential interface force t 0 in horizontal direction (b) (A), tangential interface moment m 0 in horizontal direction (B) and ratio between tangential interface force and interface wall moment C = t 0 9 d 50 /m 0 in horizontal direction versus normalized interface roughness parameter R n at grooves' height (C) at residual state (e o = 0.55 and r n = 100 kPa)

Conclusions
Comprehensive DEM simulations of the sand-structure interface behaviour with varying interface roughness were performed for cohesionless sand that was simulated by non-symmetric convex granular clumps of the same shape (aspect ratio equal to 1.5). The results were confronted with similar simulation outcomes performed for pure spheres with contact moments [15] to show the effect of the particle shape on the interface resistance and shear localization. Interface boundary conditions within micropolar continua were proposed. The following major conclusions can be summarized from our simulations: • As compared to spheres with contact moments, clumps provided in general due to a different particle shape the higher residual interface resistance, coordination numbers, particle rotations and slips, and the lower thickness of the shear zone at the interface, volume changes, orientation range of mean contact forces to the horizontal, ratios between the particle rotation multiplied by the mean grain diameter and slip along at the interface and ratios between the horizontal tangential interface force multiplied by the mean grain diameter force and horizontal tangential interface moment. The peak interface resistance solely was similar. The rolling resistance model was proved to be limited to capture particle shape effects. • The normalized interface roughness had a huge influence on the mobilized interface friction angle and thickness of the interface shear zone. The interface resistance increased with increasing interface roughness. The peak and the residual interface friction angle grew with the raising normalized roughness parameter. The relationship between the peak/residual interface friction angle and normalized interface roughness was parabolic/bi-linear. The residual interface friction angle was by 3°-7°lower than for pure spheres with contact moments (R n C 0.25). • The sand-structure interface resistance was strictly combined with the thickness of the shear zone at the interface wherein pronounced grain rotations occurred. The thickness of the shear zone at the interface increased with the growing normalized interface roughness parameter in a bi-linear way as the residual interface friction angle. The thickness of the interface shear zone expanded in a parabolic way from 1 9 d 50 up to 12 9 d 50 for R n = 0.01-2.0. It was smaller by 15-40% than for pure spheres with contact moments (R n C 0.25). • The coordination number was higher by 55-70% than for pure spheres with contact moments.
• For the very rough interface (R n C 0.75), the largest sphere rotation was located slightly above the interface (h/d 50 = 3-5) and for R n B 0.5 it was directly located at the interface where it diminished with the reduction of R n . For R n C 0.75, the sphere rotations were approaching zero at the interface (the spheres were trapped in asperities). The largest rotation at the interface was higher by 30% than for pure spheres with contact moments (R n C 0.25). • The ratio between the horizontal slip and total horizontal displacement along the interface decreased with increasing interface roughness. This ratio was similar (R n B 0.25) or slightly higher (R n [ 0.25) than for pure spheres with contact moments. The evolution was hyperbolic in the range 0 B R n B 2.0 and almost linear in the range 0 B R n B 1.0. • The ratio between the interface grain rotation multiplied by the mean grain diameter and interface grain slip was almost the same during shearing at the interface. It systematically reduced with decreasing normalized interface roughness parameter. It was twice smaller as for pure spheres with contact moments. • The evolution of the ratio between the horizontal tangential force and horizontal tangential moment at the interface changed in a hyperbolic way between 1.0 (R n = 2) and 40 (R n = 0.01) and between 1.0 (R n = 2) and 8.0 (R n = 0.10). The evolution shape was similar to this for pure spheres with contact moments. • The non-uniformity of contact forces in granular specimens at the interface increased with the growing normalized interface roughness. It was lower than for pure spheres with contact moments. The orientation of mean contact forces to the horizontal at the residual state also grew with the increasing normalized interface roughness (from 135°up to 140°). This orientation range was lower than for pure spheres with contact moments (125°-160°). • Three different relationships were proposed to describe the interface boundary conditions in micropolar continua. They included the ratio between the micropolar rotation multiplied by the mean grain diameter and slip along the interface, the ratio of the slip at the interface and total imposed horizontal displacement and the ratio between the horizontal tangential force multiplied by the mean grain diameter and horizontal tangential moment along the interface. All those ratios slightly depended on the particle shape. They were approximately related to the normalized interface roughness parameter only. 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/.