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

Abstract

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.

1 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 micro-polar 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.

2 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 micro-polar 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].

2.1 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]:

$$\vec{F}_{n} = K_{n} U\vec{N}$$

(1)

$$\vec{F}_{s} = \vec{F}_{s,prev} + \Delta \vec{F}_{s} \quad {\text{with}}\quad \Delta \vec{F}_{s} = K_{s} \Delta \vec{X}_{s} ,$$

(2)

$$K_{n} = E_{c} \frac{{2R_{A} R_{B} }}{{R_{A} + R_{B} }}\quad {\text{and}}\quad K_{s} = \upsilon_{c} E_{c} \frac{{2R_{A} R_{B} }}{{R_{A} + R_{B} }},$$

(3)

$$||\vec{F}_{s}|| -|| \vec{F}_{n}|| \times \tan \mu \le 0,$$

(4)

$$\Delta M = K_{r} \Delta \vec{\omega }\quad {\text{with}}\quad K_{r} = \beta K_{s} R_{A} R_{B} ,$$

(5)

$$||\vec{M}|| - \eta \frac{{R_{A} + R_{B} }}{2}||\vec{F}_{n} || \le 0$$

(6)

$$\vec{F}_{damp}^{k} = \vec{F}^{k} - \alpha_{d} \cdot \text{sgn} \left( {\vec{\nu }^{k} } \right)\left| {\vec{F}^{k} } \right|\quad {\text{and}}\quad \vec{M}_{damp}^{k} = \vec{M}^{k} - \alpha_{d} \cdot \text{sgn} \left( {\overrightarrow {{\dot{\omega }}}^{k} } \right)\left| {\vec{M}^{k} } \right|,$$

(7)

where $\vec{F}_{n}$—the normal contact force, U—the overlap between discrete elements, $\vec{F}_{s}$—the tangential contact force, $\vec{F}_{s,prev}$—the tangential contact force from the previous iteration, $\vec{N}$—the unit normal vector at each contact point, $\vec{X}_{s}$—the relative tangential displacement of the sphere centre, K_n—the normal contact stiffness, K_s—the tangential contact stiffness, μ—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, ν_c—the Poisson’s ratio of the grain contact, M—the contact moment, K_r—the rolling contact stiffness, β—the dimensionless rolling stiffness coefficient, $\vec{\omega }$—the resultant angular rotation between two elements, η—the dimensionless limit rolling coefficient, $\vec{F}_{damp}^{k}$ and $\vec{M}_{damp}^{k}$—the damped contact force and moment, $\vec{F}^{k}$ and $\vec{M}^{k}$—the kth components of the residual contact force and contact moment vector, $\vec{\nu }^{k}$ and $\overrightarrow {{\dot{\omega }}}^{k}$ are the kth components of the translational and rotational velocities of spheres and α_d—the 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 ν [28, 30]. The effect of damping was negligible in quasi-static calculations [28, 30].

Fig. 1

Mechanical response of linear contact model without (A) and with contact moments (A + B): a tangential contact model, b normal contact model and c rolling contact model and C) loading and unloading path (tangential and rolling contact) [15, 24, 62]

The five main local material parameters are necessary for our DEM simulations: E_c (modulus of elasticity of the grain contact), ν_c (Poisson’s ratio of the grain contact), μ (inter-particle friction angle), β (rolling stiffness coefficient) and η (limit rolling coefficient). In addition, a particle radius R, particle mass density ρ and numerical damping parameter α are required. The DEM material parameters: E_c, ν_c, μ, β, η and α were calibrated using the corresponding homogeneous axisymmetric triaxial laboratory test results on Karlsruhe sand with the different initial void ratio and lateral pressure [26, 67]. The procedure for determining the material parameters in DEM was described in detail by Kozicki et al. [28, 30]. Note that the representative elastic contact moduli E_c and ν_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 × 100 × 20 mm³. 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≤ 0.1 × d₅₀), rough (0.1 × d₅₀< R_max< 0.5 × d₅₀) and very rough (R_max≥ d₅₀), where R_max—the maximum vertical distance between peaks and valleys evaluated along the wall over the length of 3 × d₅₀ [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₅₀ = 0.50 mm, grain size between 0.08 mm and 1.8 mm, uniformity coefficient U_c = 2, maximum specific weight γ ^max_d  = 17.4 kN/m³, minimum void ratio e_min = 0.53, minimum specific weight γ ^min_d  = 14.6 kN/m³ and maximum void ratio e_max = 0.84. The sand grains were classified as sub-rounded/sub-angular.

3 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, ν_c, μ, β, η and α), 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 × 10 × 10 cm³, 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₅₀ = 5 mm (10 times larger than the real one). The mass density was 2600 kg/m³. 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 inter-particle friction angle was varied between 0° and μ (initially dense sand) and between 89° and μ (initially loose sand) to exactly reproduce the target initial void ratio. During dynamic compression to the desired confining pressure σ_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 μ. The DEM simulations were carried out for initially dense sand (initial void ratio e_o = 0.53) for three different lateral pressures σ_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].

Fig. 2

Triaxial compression test on initially dense sand (initial void ratio e_o = 0.53): A view on granular specimen between walls in DEM simulations and B vertical normal stress σ₁ and C volumetric strain ε_v versus vertical normal strain ε₁ from DEM compared to experiments [27] for different initial lateral pressures: σ_c = 50 kPa, σ_c = 200 kPa and σ_c = 500 kPa

Table 1 Material parameters assumed in all DEM simulations

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. For initially dense sand and confining pressure σ_c = 200 kPa, the global macroscopic elastic parameters were E  = 80 MPa and ν  = 0.25 using spheres with contact moments and a linear contact model, and E  = 70 MPa and ν  = 0.22 using spheres with contact moments and a nonlinear contact model. The global maximum mobilized internal friction angle φ_max and dilatancy angle were φ_max = 42^o and ψ  = 30^o (spheres with contact moments and a linear contact model), and φ_max = 42.1^o and ψ  = 30^o (spheres with contact moments and a nonlinear contact model). Due to negligible differences, further simulations were carried out with a linear contact model to reduce the computation time. It has been shown that the Hertz–Mindlin–Deresiewicz model due to its complication is computationally very time-consuming with a large number of particles in contact [41, 47, 76].

4 DEM simulations of monotonic direct wall shearing tests

4.1 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.

Fig. 3

Direct wall friction test in DEM: A geometry of three-dimensional DEM model, B sand behaviour during shearing (a) initial state and b) final state) and C rigid bottom wall sections with different normalized wall roughness parameter R_n = h_g/d₅₀: a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 (h_g—groove height and d₅₀—mean grain diameter)

The grain diameter of sand linearly varied between 0.25 mm and 0.75 mm with the mean grain diameter of d₅₀ = 0.5 mm (as in the experiment). All walls confining the sand specimen were assumed as rigid boundaries. The top and side walls were movable whereas the bottom wall was fixed. To induce shearing along the bottom (Fig. 3b), the sand was horizontally sheared under a constant velocity of the shear box (from the left to the right) until the shear displacement reached the limit of 7 mm.

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₅₀, where h_g is the groove height and d₅₀ 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 = 2×d₅₀). The groove inclination to the bottom α_g diminished with decreasing R_n (e.g. α_g = 45^o for R_n = 1.0) (Fig. 3C). Some DEM simulations were also carried with the different groove inclination to the bottom α_g and groove distance s_g for R_n = 1.0 (Sect. 5.4). The constant uniform vertical pressure σ_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⁻⁴). The wall friction angle between particles and grooves was assumed to be μ_sw = μ = 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 φ_w = arctan(T/N) versus the horizontal shear box displacement u_x for the normalized wall roughness parameter R_n = 1.0 with initial void ratio of e_o = 0.55 and pressure of σ_n = 100 kPa. The resultant normal force N constituted the sum of vertical forces acting on the top wall and the resultant friction wall force T the sum of horizontal forces acting on the top and side walls.

Fig. 4

Calculated mobilized wall friction angle φ_w versus horizontal displacement u_x from DEM for initially dense sand (initial void ratio e₀ = 0.55) with vertical pressure σ_n = 100 kPa and normalized wall roughness parameter R_n = 1.0: A) effect of 2D (a) and 3D simulations (b) and B) effect of different specimen width D: a D  = 100 mm, b D  = 50 mm, c D  = 25 mm, d D  = 5 mm and e D  = 2.5 mm

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 × d₅₀) in all DEM simulations instead of D  = 100 mm as in the experiment.

4.2 Effect of initial void ratio of sand and vertical pressure

To analyse the effect of initial void ratio e₀ on the specimen behaviour along the rough bottom wall, two different initial void ratios were assumed: e₀ = 0.55 (initially dense sand) and e₀ = 0.80 (initially loose sand) with the normalized wall roughness parameter R_n = 1.0 and vertical pressure σ_n = 100 kPa. The DEM results with e₀ = 0.55 and e₀ = 0.80 were compared with our experiments [49, 53] (Fig. 5). The DEM calculations were also carried out with three different vertical pressures σ_n: σ_n = 50 kPa, σ_n = 200 kPa and σ_n = 500 kPa (R_n = 1.0 and e_o = 0.55) (Fig. 6).

Fig. 5

DEM results of mobilized wall friction angle φ_w and volumetric strain ε_v versus horizontal displacement u_x (continuous lines) compared to experimental results [49, 53] (dashed lines) for normalized roughness parameter R_n = 1.0 with two different initial void ratios of sand: e₀ = 0.55 and e₀ = 0.80

Fig. 6

Mobilized wall friction angle φ_w (A) and volumetric strain ε_v versus horizontal displacement u_x (B) from DEM for different vertical pressure σ_n: a σ_n = 50 kPa, b σ_n = 200 kPa and c σ_n = 500 kPa (normalized wall roughness parameter R_n = 1.0 and initial void ratio of sand e_o = 0.55)

The evolutions of φ_w = f(u_x) and ε_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 φ_w,max = 49.6º at the peak (u_x = 0.5 mm) and φ_w,res = 36º (u_x> 3 mm) at the residual state for initially dense sand, and φ_w,max = 42º at the peak (u_x = 1.5 mm) and φ_w,res = 34º (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 σ_n. Globally, the initially dense sand dilated and the initially loose sand contracted (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 φ_w,max were smaller by the factor 2 than in the experiment. The calculated residual wall friction angles φ_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 σ_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 φ_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 φ_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).

4.3 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.55) and vertical pressure of σ_n = 100 kPa (Fig. 7). The values of the peak wall friction angle φ_w,max, residual wall friction angle φ_w,res and volumetric strain ε_v are given in Table 2. The relationship between the computed values of φ_w,max and φ_w,res and R_n is shown in Fig. 8. The DEM results were compared with our experiments [53] in Fig. 9.

Fig. 7

Mobilized wall friction angle φ_w (A) and volumetric strain ε_v versus horizontal displacement u_x (B) from DEM for different normalized wall roughness parameter R_n (initial void ratio e₀ = 0.55 and vertical pressure σ_n = 100 kPa): a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01

Table 2 Calculated values of peak wall friction angle φ_w,max, residual wall friction angle φ_w,res and residual volumetric strain ε_v for different normalized wall roughness parameter R_n (regular grooves)

Fig. 8

Relationship between peak wall friction angle φ_w,max and residual wall friction angle φ_w,res and normalized wall roughness parameter R_n from DEM: a) φ_w,max and b) φ_w,res

Fig. 9

DEM results of mobilized wall friction angle φ_w (A) and volumetric strain ε_v (B) versus horizontal displacement u_x (continuous lines) for initially dense specimen (e_o = 0.55) under vertical pressure of σ_n = 100 kPa compared to experimental results [49, 53] (dashed lines) for different normalized wall roughness parameter R_n: a R_n = 1.0 (very rough surface), b R_n = 0.25 (rough wall) and c R_n = 0.01 (smooth wall)

5 Evolution of mobilized wall friction angle and volumetric strain

The peak φ_w,max and the residual wall friction angle φ_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≥ 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≥ 0.75 and reduced with decreasing R_n (Fig. 7B). The relationship between the values of φ_w,max and R_n and between φ_w,res and R_n was bi-linear as, e.g. in experiments [19, 48] and DEM analyses [25, 70] (Fig. 8). The value of φ_w,max = 18º for R_n = 0.01 was obviously equal to the assumed wall friction angle between particles and grooves μ_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).

6 Distribution of particle rotation, particle displacement and void ratio

The effect of the normalized wall roughness parameter R_n on the sphere rotations ω in the entire sand specimen is shown at the residual state (u_x = 7 mm) in Fig. 10. The values of rotations ω were calculated from a cubic averaging cell of the size 5d₅₀× 5d₅₀× 5d₅₀ being moved by d₅₀. In the averaging cell, the centres of spheres were considered independently of their diameter. Figure 11 presents the zoom on single sphere rotations in the granular segment at the front side (mid-length), and Fig. 12 shows the distribution of sphere horizontal displacements u_x, sphere rotations ω and void ratio e across the normalized specimen height h/d₅₀ at the specimen mid-region at the residual state. The values of u_x, ω and e were calculated from the averaging cell of the size 5d₅₀× 5d₅₀× 1d₅₀ (length × width × height) being moved by d₅₀. 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 σ_n on the distribution of sphere rotations ω across the normalized specimen height h/d₅₀ is demonstrated in Fig. 14 at the residual state. The distribution of the ratio between the sphere rotation multiplied by the mean grain diameter ωd₅₀ and the sphere slip u [(ωd₅₀)/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 (ωd₅₀)/u against the inversed normalized wall roughness parameter 1/R_n at the grooves’ height is given in Fig. 16.

Fig. 10

Distribution of sphere rotations ω with attached scale in [rad] at residual state for u_x = 7 mm from DEM with: A different normalized wall roughness parameter R_n: a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 (e_o = 0.55 and σ_n = 100 kPa) and B different initial void ratio of sand (e_o = 0.80, σ_n = 100 kPa and R_n = 1.0) (positive sign—clockwise rotation, negative sign—anticlockwise rotation)

Fig. 11

Zoom on distribution of sphere rotations across granular segment at front side for residual state u_x = 7.0 mm from DEM: A with different normalized wall roughness parameter R_n (e_o = 0.55 and σ_n = 100 kPa): a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 and B with different initial void ratio (e_o = 0.80, σ_n = 100 kPa and R_n = 1.0) (red colour—clockwise rotations, blue colour—anticlockwise rotations, white colour—no rotation, colour scale is in [rad])

Fig. 12

Distribution of horizontal sphere displacement u (A), sphere rotation ω (B) and void ratio e (C) across normalized specimen height h/d₅₀ at the specimen mid-point at residual state for u_x = 7.0 mm from DEM with different normalized wall roughness parameter R_n: a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 (e_o = 0.55 and σ_n = 100 kPa)

Fig. 13

Relationship between thickness of wall shear zone t_s and normalized wall roughness parameter R_n (e₀ = 0.55, σ_n = 100 kPa)

Fig. 14

Distribution of sphere rotations ω across normalized specimen height h/d₅₀ at specimen mid-point at residual state for u_x = 7.0 mm from DEM (R_n = 1.0 and e_o = 0.55) with: A) different initial void ratio e₀ of sand (R_n = 1.0 and σ_n = 100 kPa): a e₀ = 0.55 and b e₀ = 0.80 and B) different vertical pressure σ_n: a σ_n = 50 kPa, b σ_n = 200 kPa and c σ_n = 500 kPa

Fig. 15

Distribution of ratio ωd₅₀/u across normalized specimen height h/d₅₀ at specimen mid-point at residual state for u_x = 7.0 mm from DEM with different normalized wall roughness parameter R_n: a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 (initial void ratio e_o = 0.55 and vertical pressure σ_n = 100 kPa)

Fig. 16

Relationship between wall grain rotation—wall grain slip ratio (ωd₅₀)/u and inversed normalized wall roughness parameter 1/R_n

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 ω where the rotation was ω ≤ 5% of the maximum sphere rotation ω_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 × d₅₀ (R_n = 0.75–2.0), t_s = 10 × d₅₀ (R_n = 0.50), t_s = 6×d₅₀ (R_n = 0.25), t_s = 2×d₅₀ (R_n = 0.10) and t_s = 1×d₅₀ (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 σ_n = 40–100 kPa [24] and for R_n = 0.02 and σ_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≥ 0.75, the largest sphere rotation was located slightly above the bottom wall (h/d₅₀ = 5–6), and for R_n≤ 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≥ 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≥ 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.55 (R_n = 0.01) and e  = 0.85 (R_n≥ 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  = ωd₅₀/u systematically reduced with decreasing R_n (Fig. 15) or increasing 1/R_n (Fig. 16): A  = 15,000 (R_n = 2.0), A  = 27.75 (R_n = 1.0), A  = 0.85 (R_n = 0.75), A  = 0.25 (R_n = 0.50), A  = 0.07 (R_n = 0.25), A  = 0.03 (R_n = 0.10) and A  = 0.01 (R_n = 0.01). The ratio ωd₅₀/u was approximately equal to R_n at the residual state [(ωd₅₀)/u ≅ R_n] and might be, thus, used as a boundary condition in micro-polar continua [55]. For the case of R_n≥ 0.75, where all sphere rotations ω tended to zero, the boundary condition could be simplified as ω  = 0 and u  = 0 (no rotations and displacements) as in the experiment [49].

7 Distribution of contact forces

Figures 17 and 18 demonstrate the distribution of the normal (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.

Fig. 17

Distribution of contact normal forces in granular specimen at residual state for u_x = 7.0 mm from DEM with: A) different normalized roughness parameter (e_o = 0.55 and σ_n = 100 kPa): a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 and B) different initial void ratio (e_o = 0.80, σ_n = 100 kPa and R_n = 1.0) (red colour corresponds to normal contact forces higher than mean value, maximum value of forces is 0.50 N)

Fig. 18

Distribution of contact shear forces in granular specimen at residual state for u_x = 7.0 mm from DEM with: A) different normalized roughness parameter (e_o = 0.55 and σ_n = 100 kPa): a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 and B) different initial void ratio (e_o = 0.80, σ_n = 100 kPa and R_n = 1.0) (red colour corresponds to normal contact forces higher than mean value, maximum value of forces is 0.15 N)

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 non-uniformity 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 σ_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).

Fig. 19

Polar mean contact force distribution in granular specimen at beginning of test (black line) and at residual state (red line) from DEM with: A) different normalized roughness parameter (e_o = 0.55 and σ_n = 100 kPa): a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 and B) different initial void ratio (e_o = 0.80, σ_n = 100 kPa and R_n = 1.0)

The number of all sphere contacts, expressed by the coordination number c (Fig. 20), was correlated with the change of volumetric strain ε_v and void ratio e (Figs. 5B and 6B). For initially dense sand (R_n = 1.0 and σ_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.

Fig. 20

Evolution of coordination number c in initially dense sand (e₀ = 0.55) from DEM simulations (normalized wall roughness parameter R_n = 1.0 and vertical pressure σ_n = 100 kPa)

8 Distribution of stresses and couple stresses

The stresses and couple stresses of a single sphere were calculated as [34]

$$\sigma_{ij} = \frac{1}{{V^{p} }}\mathop \sum \limits_{c = 1}^{N} x_{i}^{c} f_{j}^{c} ,$$

(8)

and

$$m_{ij} = \frac{1}{{V^{p} }}\mathop \sum \limits_{c = 1}^{N} x_{i}^{c} m_{j}^{c} ,$$

(9)

where N—the contact points, $x_{i}^{c}$—the ith component of branch vector jointing from the centre of mass of the particle to the contact point ‘c’, $f_{j}^{c}$—the jth component of the total force at the contact point ‘c’, $m_{j}^{c}$—the jth component of the total moment at the contact point ‘c’ and V^p—the cell volume. The values of σ_ij and m_ij were again calculated from an averaging cell of the size 5d₅₀× 5d₅₀× 1d₅₀ being moved by d₅₀.

Figures 21 and 22 show the distributions of the normal stresses σ_ii, shear stresses σ_ij and stress moments m_ij across the normalized specimen height h/d₅₀ at the specimen mid-point at residual state for u_x = 7.0 mm for the different normalized wall roughness parameter R_n (σ₁₁—horizontal normal stress, σ₂₂—vertical normal stress, σ₁₂—horizontal shear stress, σ₂₁—vertical shear stress, m₂₃—horizontal stress moment and m₁₃—vertical stress moment). The relationship between the wall shear stress—wall couple moment ratio (σ₁₂× d₅₀)/m₂₃ and the normalized wall roughness parameter R_n [55] at the wall in the residual state is demonstrated in Fig. 23.

Fig. 21

Distribution of horizontal σ₁₂ (A) and vertical shear stress σ₂₁ (B) and vertical m₁₃ (C) and horizontal stress moment m₂₃ (D) across normalized height h/d₅₀ at residual state at specimen mid-point with different normalized wall roughness parameter R_n: a R_n = 2.0, b R_n = 1.0, c R_n = 0.75, d R_n = 0.50, e R_n = 0.25, f R_n = 0.10 and g R_n = 0.01 (e_o = 0.55 and σ_n = 100 kPa)

Fig. 22

Distribution of stress components: horizontal normal stress σ₁₁ (a), vertical normal stress σ₂₂ (b), horizontal shear stress σ₁₂ (c) and vertical shear stress σ₂₁ (d) across normalized specimen height h/d₅₀ at residual state at specimen mid-point with different normalized wall roughness parameter R_n: A) R_n = 1.0, B) R_n = 0.50 and C) R_n = 0.01 (e_o = 0.55 and σ_n = 100 kPa)

Fig. 23

Relationship between wall shear stress— wall horizontal couple stress ratio (σ₁₂× d₅₀)/m₂₃ and normalized wall roughness parameter R_n in residual state (e_o = 0.55 and σ_n = 100 kPa)

The wall shear stresses σ₁₂ and σ₂₁ grew with increasing wall roughness and were non-symmetric at the wall (σ_12≠σ₁₂) (Figs. 21 and 22). The distribution of the vertical normal stress σ₂₂ was almost the same across the specimen height in contrast to the horizontal normal stress σ₁₁ and shear stresses σ₁₂ and σ₂₁. The wall couple stresses m₁₃ (vertical) and m₂₃ (horizontal) diminished with reducing wall roughness in the range 0.25 ≤ R_n≤ 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 (σ₁₂× d₅₀)/m₂₃ and R_n was solely obtained for 0.25 ≤ R_n≤ 2.0 (Fig. 23).

9 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 mid-region at the residual state for u_x = 7.0 mm with the different normalized wall roughness parameter R_n (e_o = 0.55 and σ_n = 100 kPa). The wall values were computed from the area 5d₅₀× 5d₅₀. 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’ × d₅₀)/m’ in the global horizontal direction acting at the groove height on the granular segment.

Fig. 24

Normal and tangential wall forces n’ and t’ and tangential wall moments m’ acting on inclined wall grooves

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’ × d₅₀)/m’ in horizontal direction versus normalized roughness parameter R_n at grooves’ height (C) at residual state (e_o = 0.55 and σ_n = 100 kPa)

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 bi-linear (Fig. 25Ab) as in the case of the wall friction angles φ_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’ × d₅₀)/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’ × d₅₀)/m′  = const might be used as a wall boundary condition at wall nodes in FEM using micro-polar continua.

9.1 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 × d₅₀ and 16 × d₅₀) and grooves’ inclination to the bottom α_g (45^o or 27^o) (Fig. 3) were analysed in detail (Figs. 26, 27 and 28).

Fig. 26

Various types of wall grooves assumed in DEM simulations (R_n = 1.0): A type ‘1’ (basic with (h_g = d₅₀, α_g = 45^o and s_g = 2×d₅₀), B type ‘2’ (h_g = d₅₀, α_g = 45^o, s_g = 4×d₅₀), C type ‘3’ (h_g = d₅₀, α_g = 45^o, s_g = 8×d₅₀), D type ‘4’ (h_g = d₅₀, α_g = 45^o, s_g = 16 × d₅₀) and E type ‘5’ (h_g = d₅₀, α_g = 27^o, s_g  = 4×d₅₀)

Fig. 27

Mobilized wall friction angle φ_w (A) and volumetric strain ε_v versus horizontal displacement u_x (B) from DEM for different wall groove types of Fig. 23 (R_n = 1.0, σ_n = 100 kPa and e_o = 0.55): a type ‘1’, b type ‘2’, c type ‘3’, d type ‘4’ and e type ‘5’

Fig. 28

Distribution of sphere rotations ω in entire granular specimen with attached scale in [rad] at residual state for u_x = 7 mm for various types of wall grooves of Fig. 23 from DEM (e_o = 0.55 and σ_n = 100 kPa): A type ‘1’ (h_g = d₅₀, α_g = 45^o, s_g = 2×d₅₀), B type ‘2’ (h_g = d₅₀, α_g = 45^o, s_g = 4×d₅₀), C type ‘3’ (h_g = d₅₀, α_g = 45^o, s_g = 8×d₅₀), D type ‘4’ (h_g = d₅₀, α_g = 45^o, s_g = 16 × d₅₀) and E type ‘5’ (h_g = d₅₀, α_g = 27⁰, s_g  = 4×d₅₀)

The peak wall friction angle, maximum volumetric strain and thickness of the wall shear zone decreased with increasing s_g and with decreasing α_g (Figs. 27 and 28). The effect of the groove distance s_g proved to be small and the effect of the groove inclination α_g was strong on the peak wall friction angle (Fig. 27). The effect of the groove distance started to be visible for s_g≥ 8×d₅₀. The peak wall friction angle was lower by 5% for s_g≥ 8×d₅₀ than for s_g = (2–4) × d₅₀. The asperity angle of α_g = 45^o yielded the higher maximum wall friction angle than α_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].

10 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 × d₅₀ up to 14 × d₅₀ 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 α_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.

—The horizontal grain slip along the bottom obviously decreased with increasing normalized wall roughness parameter R_n 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 at the residual state.

—For very rough walls (R_n≥ 0.75), the largest sphere rotation was located slightly above the bottom wall (h/d₅₀ = 5–6), and for R_n ≤ 0.5, it was directly located at the bottom wall where it diminished with the reduction of R_n. For R_n ≥ 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 ωd₅₀/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’ × d₅₀)/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 initially dense sand) and diminishing initial void ratio (from 115^o up to 160^o for a very rough wall).