Effect of grain roughness on strength, volume changes, elastic and dissipated energies during quasi-static homogeneous triaxial compression using DEM

A quasi-static homogeneous drained triaxial compression test on cohesionless sand under constant lateral pressure was simulated using a three-dimensional DEM model. Grain roughness was modelled by means of symmetric clusters composed of rigid spheres imitating irregular particle shapes. The effect of grain roughness on shear strength, dilatancy, kinetic, elastic and dissipated energies was numerically analyzed. Some numerical results were compared with available experimental results.

and hypoplastic constitutive models enhanced by a characteristic length of micro-structure to describe strain localization (e.g. [7,13,[38][39][40][41]). In turn, the latter ones perform simulations at the grain scale, i.e. each grain is modelled individually (e.g. [2,4,16,24,30,31,42]). Their advantages are that they directly model micro-structure and can be used to comprehensively study the mechanism of the initiation, growth and formation of shear zones at the micro-level which strongly affect macro-properties of granular matter. The disadvantages are: high computational cost, inability to model grain shape accurately and difficulty to validate it experimentally as the inertial and damping effects lose their meaning in quasi-static problems. However, they become more and more popular nowadays for modelling granular materials due to an increasing speed of computers and a connection possibility to the finite element method [33].
The concept of the total and stored elastic energy in soil mechanics was already discussed in the literature (see [8] and references therein). In the case of metals the elastic energy is divided into fully recoverable and stored or hidden portions. During unloading the recoverable part is released but the hidden energy remains stored in the material in the unloaded state as it represents the energy of residual microstresses at the grain and crystalline lattice levels. In soils and especially in granular materials the process of full unloading releases totally the grain interaction and related elastic energy. However, the partial unloading with the remained hydrostatic pressure can be considered and then the division into recoverable and stored elastic energies can be assumed. The free energy of a representative soil element can then be decomposed into elastic strain and irreversible (plastic) strain dependent portions = e (ε e ) + p (ε p ) [9]. Here the first term represents the elastic strain dependent energy for a granular aggregate when only elastic grain deformation occurs with no sliding at contact interfaces and the second term represents the stored energy due to contact slip and sliding. It is not clear whether the concept of separation of energies is valid for granular materials since it is difficult to execute an elastic loading or unloading process for which the contact force interaction would occur only along normal directions to contact interfaces. In our study the evolution of the total elastic energy will be analyzed during the continuing deformation process without separation into recoverable and stored terms.
The objective of this paper is to present numerical analyses of quasi-static homogeneous true triaxial compression tests carried out to determine the macroscopic behaviour of sand specimens composed of discrete elements in the form of clusters. A three-dimensional discrete model YADE developed at University of Grenoble was used [21,22]. The particle breakage was not considered. The discrete simulation results were compared with the corresponding experimental data from drained axisymmetric triaxial compression tests performed by [45] at Karlsruhe University with real sand. The intention of our studies was to calculate the effect of the grain roughness (shape) on the shear strength, dilatancy, elastic and dissipated energies of real sand (so-called Karlsruhe sand), which had the same initial void ratio, mean grain diameter and grain distribution. A special attention was paid to the energy transformation in sand and its elastic and dissipative characteristics, playing a fundamental role in the granular matter behaviour [1,2,37,46]. The energy and dissipation results were compared with the similar ones from simplified two-dimensional simulations of biaxial compression with round particles performed by [23] and by [6]. In addition, the original own method was introduced to generate different grain shapes in the form of clusters of spheres.

Discrete model
The discrete element method (DEM) is widely used to model a range of processes across many industries [2,10,14,19,23,32,42,47]. To simulate the behaviour of sand, a threedimensional spherical discrete model YADE was developed at University of Grenoble [21,22] by taking advantage of the so-called soft-particle approach (i.e. the model allows for particle deformation which is modelled as an overlap of particles). A dynamic behaviour of the discrete Fig. 1 Tangential and normal contact model [21,36] system is solved numerically using a force-displacement Lagrangian approach and tracks the positions, velocities, and accelerations of each particle individually. It uses an explicit finite difference algorithm assuming that velocities and accelerations are constant in each time step. To calculate forces acting in particle-particle or particle-wall contacts, a particle interaction model is assumed in which the forces are typically subdivided into normal and tangential components. The total forces acting on each particle are summed. Next, the problem is reduced to the integration of Newton's equations of motion for both translational and rotational degrees of freedom. As the results, accelerations of each particle are obtained. The time step is incremented and accelerations are integrated over time to determine updated particle velocities and positions. To maintain the numerical stability of the method and to obtain a quick convergence to a quasi-static state of equilibrium of the assembly of particles, damping forces are introduced [10].
The interaction force vector F representing the action between two spherical discrete elements in contact is decomposed into a normal and tangential vector, respectively. A linear elasticity is chosen in our contact model. The normal and tangential forces are linked to the displacements through the normal stiffness K n and the tangential stiffness K s (Fig. 1a, b) where U is the penetration depth between elements, → N denotes the normal vector at the contact point and X s is the incremental tangential displacement. The tangential force F s is obtained by summing its increments. The stiffness parameters are calculated with the aid of the modulus of elasticity of the grain contact E c and two neighbouring grain radii R A and R B (to determine the normal stiffness K n ) and with the aid of the modulus of elasticity E c and Poisson's ratio ν c of the grain contact and two neighbouring grain radii R A and R B (to determine the tangential stiffness K s ), respectively [36] Fig. 2 Experimental results of axisymmetric triaxial tests by [45] with Karlsruhe sand (mean grain diameter d 50 = 0.5 mm): relationship between σ 1 /σ c and ε 1 , and between ε v and ε 1 at different confining pressures σ c : a initially very dense (e o = 0.53) sand, b initially loose sand (e o = 0.80) (σ 1 -vertical normal stress, If the grain radius R A = R B = R, the stiffness parameters are equal to: K n = E c R and K s = ν c E c R (thus K s /K n = ν c ), respectively. The frictional sliding starts at the contact point if the contact forces F s and F n satisfy the frictional Mohr-Coulomb equation (Fig. 1a) with μ as the inter-particle friction angle (tension is not allowed). No forces are transmitted when grains are separated. The assumed tangential and normal contact relationships in the discrete model are demonstrated in Fig. 1a, b, respectively (the unloading is purely elastic). The elastic parameters K n and K s play a major role only in the elastic response of granulates. The normal stiffness modulus K n is related to the average modulus of elasticity of granular material E and the Poisson ratio ν. The ratio K s /K n = ν c depends on the Poisson's ratio ν (the relationship is approximately linear, i.e. E R ≈ K n /5, ν ≈ ν c /2, cf. [5]). A choice of a linear elastic normal contact ( Fig. 1b) is not in agreement with a non-linear Hertz interaction law between two adjacent elastic spheres with F n ∼ U 3/2 which is the exact elastic solution [15,27]. Therefore the elastic constants of the grain contact in our model with irregularly-shaped grains do not correspond to the elastic constants of the spheres' material in the Hertz law (thus K n in Eq. 3 is several times larger than the mean normal stiffness of the spherical grain material). The elastic contact moduli are specified from the experimental data of a triaxial compression sand test, as described later in the text.
To dissipate the excessive kinetic energy in the discrete system, a simple local non-viscous damping scheme was adopted, proposed by [11], which assumes a decrease of forces which increase particle velocities by using the damping parameter α where F k is the kth component of the residual force vector and v k is the kth component of the translational velocity [36]. A positive damping coefficient α is smaller than 1 (sgn (•) returns the sign of the kth component of velocity). The equations are separately applied to each k-th component of a 3D vector x, y and z. Note that the effect of damping is insignificant in quasi-static calculations.
The following 3 main local material parameters are needed for discrete simulations:E c , ν c and μ. In addition, the particle radius R, particle density ρ and damping parameter α are required. The material parameters can be calibrated with corresponding axisymmetric triaxial laboratory test results on Karlsruhe sand by [20] and [45], In numerical simulations, a cubical sand specimen of 10× 10 × 10 cm 3 was used. A simplified linear grain distribution curve was used for Karlsruhe sand (grain range among 2.5-7.5 mm, Fig. 3). To save the computation time, discrete simulations showing the capabilities of DEM were carried out with d 50 = 5 mm (Fig. 3) instead of d 50 = 0.5 mm. The test was modelled using confining smooth rigid wall elements (without inducing shear localization). The top and bottom boundaries moved vertically as loading platens under straincontrolled conditions to simulate the confining pressure p. The loading speed was slow enough (10 mm/s) to ensure the test was conducted under quasi-static conditions. Different methods can be used to model irregularlyshaped grains [3,12,17,25,26]. In our method, each granular assembly was prepared by putting clusters of spheres  Fig. 3 (without gravity) into a cubical container with 6 external walls which had a regular cubical grid with a particle distance of 10 mm. The initial configuration of the specimen was thus isotropic. In order to obtain a desired initial density owing to grain overlapping, the inter-particle friction coefficient μ was varied between 0 and 30 • to exactly reproduce the target initial sand volumetric weight. 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 random. The assembly was then allowed to settle to a state where the kinetic energy was negligible and the friction coefficient was set to μ = 30 • . The isotropic assembly was then subjected to the boundary driven triaxial compression. In general, arbitrary symmetric and non-symmetric shapes can be obtained with the different aspect, angularity and convexity index. Figure 4 presents 12 different symmetric clusters of spheres used in discrete calculations. To generate e.g. clusters of 2 spheres (shapes 'b-f' of Fig. 4), the distance between spheres was set to be between 0 and −0.33d. In turn, the clusters of 4 and 6 spheres (shapes 'g, h' of Fig. 4) were made from 4 or 6 symmetrically placed spheres being tangential to the central point of the cluster. The discs (shape 'l' of Fig. 4) were created from 12 clumps of the type 'i' rotated around the main disc axis. Some shape indexes for the grain clusters of Fig. 4 are given in Table 1. The aspect index was defined as the ratio between the maximum and minimum cluster diameter, the convexity index '1' as the ratio between the smallest sphere volume encompassing the cluster and the cluster volume, and the convexity index '2' as the ratio between the smallest convex volume encompassing the cluster and the cluster volume.
In the case of the cluster of 2 spheres without the overlap (shape 'f' of Fig. 4), 26,300 clusters were composed of 52,600 spheres. In turn, 28,250 clusters were used with 197,750 spheres to model simple ellipsoids (shape 'i' of Fig. 4) and 14,500 clusters were used with 594,500 spheres to model disks (shape 'l' of Fig. 4). The computation time was 1 day (spheres), 2 days (cluster of 2 spheres), 4 days (ellipsoids) and 10 days (discs) using PC 3 GHz.

Discrete results of homogeneous triaxial compression test
The following discrete material parameters were used in simulations: Similarly as in the real experiment (Fig. 2), the initially dense specimens exhibits initially elasticity, hardening (connected first to contractancy and then dilatancy), reaches a peak at about of ε 1 = 2−3.5 %, gradually softens and dilates reaching at large vertical strain of 25-30 % the same value of the vertical normal stress with the specimen deforming at constant volume, i.e. a critical state is always reached. Thus, the particle shape is not essential for the global critical internal friction angle (except of the case with single spheres). The both mobilized strength and dilatancy increase in general with increasing grain roughness and rolling resistance combined with an increase of the sphere number. Thus, the irregularly shaped particles provide obviously higher internal friction angles and have less tendency to rotate than perfect circular particles. The global maximum mobilized internal friction angle increases from φ max = 28 • (spheres) up to φ max = 48.9 • (ellipsoids 'k' of Fig. 4), respectively (Fig. 6). In turn, the global residual internal friction angle increases from φ cr = 25 • (spheres 'a' of Fig. 4) up to φ cr = 32 • (disks 'l' and ellipsoids 'j' of Fig. 4), respectively (Fig. 5). The material dilatancy (volume increase) is the smallest with spheres and the highest with clumps 'f-k' of Fig. 4, respectively. The global elastic modulus is similar independently of the grain roughness (E = 70 MPa with ν = 0.25). The granular system shows small fluctuations in the residual phase. In general, the increase of the maximum internal friction angle is not directly connected with increasing three shape indexes above 1 of Table 1. Thus, the determination of the key geometrical grain parameter controlling the mechanical response (strength and volume changes) requires further investigations.    (Fig. 2). The softening rate is similar. However, the calculated strain ε 1 = 3 % is about twice too small. It might be improved by decreasing the modulus of elasticity of the grain contact [5,43] or by introducing a non-linear contact model by Hertz-Mindlin-Deresiewicz (based on our preliminary calculations). The calculated dilatancy angles of ψ = 30 • − 40 • match well the experimental outcome of ψ = 28.5 • (Fig. 2). The calculated global residual internal friction angle is 30 • − 34 • . The numerical effect of initial void ratio e o and lateral pressure σ c on the sand behaviour using clusters of 6 spheres 'h' is described in Figs. 7 and 8

Effect of grain roughness on elastic and dissipated energies
In a granular system there exist 3 main energies: elastic energy, kinetic energy and dissipated energy. In addition, the numerical dissipation also exists (see Eq. 5). The elastic internal energy stored at contacts between grains E e is expressed in terms of work of elastic contact tangential forces F s on tangential elastic displacements U t and of contact normal forces F n on penetration depths U . In general, the elastic internal energy is expressed as follows (here N -the contact number) The kinetic energy E c of grains is caused by their translation and rotation where m is the mass and I denotes the moment of inertia of a particle (v-transitional velocity, • ω-rotational velocity).
Due to quasi-static conditions, the effect of the kinetic energy E c is negligible (smaller than 1 % of the elastic energy). The dissipated energy D p is expressed in terms of work of tangential (shear) forces on related sliding displacements (see Fig. 1a) In addition, the numerical dissipation D n is specified during translation (see Eq. 5).
The total accumulated energy E = E e + E c + D p + D n is equal to the external boundary work W done on the assembly by 6 external forces on displacements of 6 rigid external walls  Fig. 4) were compared. In turn, the evolution of the external energy derivative δ E, elastic internal energy derivative δ E e and plastic dissipation derivative δ D p with respect to the vertical normal strain ε 1 is demonstrated in Fig. 10. Figure 11 shows the initial evolution of the elastic internal energy derivative δ E e with respect to the vertical normal strain ε 1 (δ E e /δε 1 ) for clusters of 2 spheres. Finally, Fig. 12 demonstrates the evolution of the kinetic energy of granular systems. There exists a roughly linear relationship between the total energy and plastic damping against the vertical normal strain (Fig. 9). The plastic dissipation during frictional sliding is equal to 50 % of the total energy at ε 1 = 3 % (corresponding to the maximum vertical stress). At the residual state of ε 1 = 30 %, it is already equal to 88 % of the total energy. The numerical damping is equal always to 6 % of the total energy. The evolution of the total elastic internal energy (and its components) is similar to the evolution of the shear strength; however the maximum value is at ε 1 = 5% (Fig. 6). At the beginning of deformation at ε 1 < 1 % (when the specimen is in the elastic range), the total energy is almost fully converted into the elastic energy. The derivative δ E e of the elastic internal work with respect to ε 1 (δ E e /δε 1 ) is initially positive (Fig. 10). It rapidly approaches zero, then becomes negative (since E e = f (ε 1 ) diminishes, Fig. 9II) at about ε 1 = 5 % (Fig. 10) and afterwards slightly increases approaching an asymptote at zero (since E e = f (ε 1 ) reaches a residual state, Fig. 9II). Beyond strains of ε 1 = 5 %, almost the entire input work is dissipated due to plastic deformation and numerical damping (the external energy rate and dissipation rate are equal δW ∼ = δ D). The energies and fluctuations of the energy rates increase with increasing grain roughness. The elastic energy decreases during a dilative deformation process only at ε 1 > 1 % (Fig. 11).
The elastic internal work is 80 % at ε 1 = 1 %, 40 % at ε 1 = 3 and 5 % at ε 1 = 30% of the total energy, respectively. The residual elastic internal work is performed by contact normal forces in 70 % and contact tangential forces in 30 %. Thus, the largest internal work is performed by the contact normal forces (due to the lack of the plastic deformation). The elastic energy ratio is the same at the residual state.
The evolution curves in Fig. 9 are qualitatively similar to those demonstrated by [6]. In turn, the evolution curves in Fig. 10 are slightly different in the initial phase only than those shown by [23], who used periodic boundary conditions (instead of walls). The calculated energy quantities are different than in analyses by [6] using the software of PFC2D, wherein e.g. the calculated elastic energy was significantly higher: 90 % (at ε 1 = 3 %) and 20 % (at ε 1 = 5 %) of the total energy (probably due to the lack of the third dimension).
The kinetic energy is very small due to the quasistatic loading of the granular system (Fig. 12). A release of the elastic energy drives grains to move. At the elastic stage, the kinetic energy is close to zero. The translational kinetic energy increases up to ε 1 =5 − 7 % and then slightly decreases. The rotational kinetic energy continuously increases (slightly). The kinetic energy show fluctuations at the residual phase (2 spheres) or already after the peak (6 spheres) which correspond to the evolution of the elastic energy and damping rate. The fluctuations increase with the grain roughness. The transitional kinetic energy is 2-5 times higher than the rotational one.

Conclusions and future work
The numerical simulations of a homogeneous triaxial compression test show that a discrete model is capable to reproduce the most important macroscopic properties of cohesionless granular materials without it being necessary to describe the granular structure perfectly. Comparing the numerical simulations with the experimental axisymmetric triaxial tests conducted for different initial void ratios and confining pressures shows that the discrete model is able to realistically predict the experimental results for real cohesionless sand.
The following detailed conclusions can be also drawn: • The model is capable of closely reproducing the behaviour of cohesionless soils in the elastic, contraction and dilatancy phase and at the critical state. At large strains, the granular specimen reaches always a critical state independently of its initial density. The higher the confining pressure, the smaller are both the global friction and dilatancy. • The sand grain roughness can be modelled by irregularlyshaped grains which may cause an increase of the strength and volume changes. • At the elastic stage, the boundary external work is mainly converted into the elastic energy. At the residual state, the almost total external boundary work is dissipated by plastic deformation. • The evolution of the elastic energy is inherently related to the dilation effect. In fact, the dilation reduces the normal contact forces and also the number of contact points. As the major contribution to the elastic energy is related to the normal forces, the elastic energy decreases during the dilative deformation process and tends to a steady state corresponding to a critical state condition. • The kinetic energy shows fluctuations which correspond to the evolution of the elastic energy and damping rate. The fluctuations increase with the grain roughness and can appear already in the softening phase if the grains are rough enough. The transitional kinetic energy is higher than the rotational one.
Our research will be continued. The discrete simulations will be carried with sand during biaxial compression out by taking into account shear localization. The local phenomena occurring in a shear zone (such as buckling of granular columns, vortices, force chain cycles, periodic alternating dilatancy and contractancy) will be carefully studied.