Derivation of Cyclic Stiffness and Strength Degradation Curves of Sands through Discrete Element Modelling

: Cyclic degradation in fully saturated sands is a liquefaction phenomenon characterized by the progressive variation of the soil strength and stiffness that occurs when the soil is subjected to cyclic loading in undrained conditions. An evaluation of the relationships between the degrada ‐ tion of the soil properties and the number of loading cycles is essential for deriving advanced cyclic constitutive soil models. Generally, the calibration of cyclic damage models can be performed through controlled laboratory tests, such as cyclic triaxial testing. However, the undrained response of soils is dependent on several factors, such as the fabric, sample preparation, initial density, initial stress state, and stress path during loading; hence, a large number of tests would be required. On the other hand, the Discrete Element Method offers an interesting approach to simulating the com ‐ plex behavior of an assembly of particles, which can be used to perform simulations of geotechnical laboratory testing. In this paper, numerical triaxial analyses of sands with different consistencies, loose and medium ‐ dense states, were performed. First, static triaxial testing was performed to char ‐ acterize the sand properties and validate the results with the literature data. Then, cyclic undrained triaxial testing was performed to investigate the impact of the number of cycles on the cyclic degra ‐ dation of the soil stiffness and strength. Laws that can be used in damage soil models were derived.


Introduction
Saturated cohesionless soils, when subjected to rapid cyclic loading in undrained conditions, manifest a progressive variation of their stiffness and strength. This phenomenon is caused by the reorganization of the complex assembly of grains which tend to compress or dilate according to their initial state of compaction. If the tendency is to compress, the pore water, that cannot flow out, will contrast this change by increasing its pressure. Consequently, the contact forces between the grains will be reduced, causing soil degradation. If the excess pore pressure is high enough to reduce each contact between the particles completely, soil liquefaction will occur. Even without reaching full liquefaction, the progressive degradation of the soil properties, during the cyclic loading, strongly affects the dynamic response of their structures [1]. Therefore, an assessment of the strength and stiffness degradation of sands is essential for the calibration of cyclic damage models [2]. The simulation of this fundamental phenomenon has been extensively studied over the last decades, and a great variety of continuum constitutive models describing nonlinear aspects of soils has been proposed [3]. Because of the continuum constitutive nature of these models, an empirical or phenomenological calibration of the model parameters is required: typically, experimental triaxial tests are performed [4][5][6].
On the other hand, the Discrete Element Method (DEM) [7] was originally developed for soil mechanics to simulate granular assemblies with individual grains modeled by disks or spherical particles [8]. Therefore, the advantage of the DEM is the ability to simulate the material microstructure and to represent directly and intrinsically the heterogeneous discrete nature of granular soils (e.g., particle shape and size distribution) [9].
DEM simulations have been conducted to replicate geotechnical testing such as the conventional monotonic triaxial test to derive the stress-strain behavior of the numerical specimen [10] and to show typical theoretical aspects such as dilatation and critical state behavior [11][12][13][14][15][16][17][18][19][20][21][22]. Most studies are based on the use of simplified sphere particle shapes to keep calculation costs as low as possible [11,14,[17][18][19][20][21][22]. On the other hand, spherical shapes can manifest excessive rolling, leading to an underestimation of the value of the friction angle as compared to real cohesionless soil [20][21][22]. To improve the reliability of the DEM model, an additional contribution, i.e., the rolling stiffness with an elastic contact law, has been introduced by Iwashita and Oda [21,22]; the resultant contact rolling moment acting against the relative rolling rotation of the particles increases the soil strength. This rolling resistance model has been used in [10] to simulate triaxial testing and the results evidenced the ability of the numerical model to reproduce complex behaviors, as observed in real granular materials, as the non-associative flow rule [23]. Moreover, rolling resistance can be an alternative solution to consider the particles' shapes instead of using computationally expensive polygons and polyhedrons [20,24]. Several studies have been based on the investigation of the irregularity of the sand particles' shape on the macroscopic behavior of the granular material [8][9][10]13,24]. In [24], a coefficient of the rolling friction related to the normalized average eccentricity of contact has been introduced to capture the effects of the grains' shape. Aboul Hosn et al. [20] performed an extensive investigation of the influence of the contact rolling resistance properties on the macroscopic behavior by performing numerical triaxial tests of sands. They showed that increasing the rolling elastic stiffness will cause an increase in the macroscopic internal friction angle; if the rolling stiffness coefficient is high enough, its influence on both the peak friction angle and the dilatancy angle becomes negligible.
The most widely used laboratory procedure to evaluate the liquefaction characteristics of sands is the cyclic triaxial test on laboratory-prepared samples [25]. The DEM has also been used to simulate fully saturated soil in undrained conditions by adopting the constant volume method [26][27][28][29][30][31]. The representation of undrained conditions through keeping volume change conditions has the advantage of saving computation time by eliminating complex fluid-particle coupling calculations. However, almost all DEM simulations of liquefiable soils focused on the micromechanical investigation of liquefaction triggering [27][28][29][30][31]. Martin et al. [28] focused on the micromechanical investigation of liquefaction performed by DEM simulations of cyclic undrained triaxial tests. Vinod et al. [31] claimed that a unique relationship exists between shear strain and the number of cycles for triggering the initial liquefaction.
Although the potentiality of the DEM simulation in capturing the behavior of real soils has been largely verified, the results have been rarely used to calibrate soil constitutive models to be used in practical applications, usually conducted through the Finite Element Method. Therefore, this paper aims to propose, for the first time, a derivation of the strength and stiffness degradation relationships to be used in cyclic damage models at a meso-or macro-scale, by exploiting the DEM instead of carrying out conventional but cumbersome and costly experimental testing. Numerical cyclic triaxial testing is, hence, performed, and the maximum deviatoric stress, as well as the shear modulus, is monitored with the increase of the number of loading cycles to derive these fundamental and practical relationships.

Discrete Element Method for Cyclic Triaxial Testing
The Discrete (or Distinct) Element Method pioneered by Cundall [7] is based on the numerical modelling of the contact interactions of an assembly of discrete particles that move under specifically prescribed boundary conditions. Essentially, the method entails computing the interaction forces generated by an artificial interpenetration of the particles and, subsequently, using these forces to calculate the new position of each particle through Newton's second law. A summary of the method is presented below.

Equation of Motion
Let us consider an assembly of particles, hereinafter embodied by circular or spherical elements, as in Figure 1. The translational motion of the i-th discrete element is governed by the Newton-Euler rigid body dynamics equations as follows: where is the acceleration of the center of mass of the particle, is its mass, and is the resultant force vector expressed as follows: , In Equation (2), is the vector of the contact forces, acting at the set of the contact points, 1 , 2 … , with each one corresponding to the middle of the artificial overall zone between any two particles (see Figure 1); is the external load vector, acting on the i-th particle at the time t, and is the global damping force vector. It is worth mentioning that represents the body forces, such as the gravitational force component • , where g is the gravity acceleration; this contribution can be neglected for quasi-static triaxial analysis, and they will not be considered in this study. The contact force vector at the contact point, in Equation (2), can be decomposed into the normal and tangential or shear components, and , respectively: The tangential or shear force, , contributes to the sliding motion between two interacting particles, as well as to the rotation about its center of mass; this rotation allows the rolling motion of the particles relative to each other, which is a fundamental aspect for the correct simulation of the redistribution of the grains in a soil sample [12,19,20]. Therefore, the rotational motion of a spherical particle of radius can be described by the following expression: in which is the angular velocity of the center of mass of the particle, is the moment of inertia tensor of components 2 5 • • 2 about its center of mass, and is the resultant moment. The resultant moment, in Equation (4), is derived from the cross product ( ) of the tangential contact forces, , by the unit normal vector, , perpendicular to the contact's tangent plane, as follows: where is the vector connecting the center of the i-th particle to the contact point , , is the rolling moment vector at the same contact points, and is the global damping moment.
Equations (1) and (4) can be solved through an explicit time integration scheme, such as the "leap-frog" algorithm [32], obtained by modifying the Verlet difference scheme [33]. Therefore, at each time, , new particle positions are computed, and a new contact or collisions are detected or updated; this generates new interaction forces that will be used in Equations (1)-(5) to reach a new state of equilibrium. In this paper, the simulation of the cyclic triaxial testing was performed under quasistatic conditions; therefore, to facilitate a rapid convergence, a purely numerical global non-viscous damping was used [34], where is the positive numerical damping coefficient, and sgn(•) returns the sign of the i-th component of the translational velocity, , and angular velocity, .

Contact Model
The Cundall-Strack elastic perfectly brittle contact model [34] was used to determine the normal and tangential contact forces in Equation (3). Therefore, the normal force acting at the i-th particle induced by the artificial overlap with the j-th particle, , is given by the following expression: (8) where is the normal stiffness and is a unit normal vector that is perpendicular to the contact's tangent plane. The shear force vector, obtained by summing up all the incremental tangential contributions, is given by the following expression: where is the shear stiffness related to the normal stiffness through the Poisson's ratio, , and Δ is the relative tangential displacement at the contact point. The incremental formulation is required to implement the Mohr-Coulomb rupture criterion, which characterizes the typical behaviour of non-cohesive geomaterials: where is the interparticle angle of friction, and ‖•‖ is the norm operator. The rolling moment, , , in Equation (5) represents the rolling resistance against the relative rolling rotation, , of two particles; it can be computed as the sum of the incremental contributions expressed as follows: in which is the rolling stiffness defined as follows: In Equation (12), is the radius of the i-th particle, and is the rolling stiffness coefficient. As for the shear force, the rolling resistance can be expressed through the frictional law: where is the limiting rolling coefficient.

Calibration of the Input Parameters for Cyclic Triaxial Testing
The contact model described in Section 2.2 can be fully determined by six parameters: the mass of each particle, , the normal stiffness, , the Poisson's ratio, , the interparticle angle of friction, , the rolling stiffness coefficient, , and the limiting rolling coefficient, . It should be noted that these parameters (except for the particle mass) do not have a direct physical relation with the microscopic properties of a real soil grain. Nevertheless, as observed in several past studies [12][13][14][15]20], a few of these parameters can be considered to have a secondary impact on the macroscopic behavior of an assembly of particles [12,20] if a few dimensionless indices are opportunely calibrated; for instance, it can be defined a dimensionless stiffness level [20]: where is the confining pressure, which in our case, should be set greater than 1000 to have a negligible effect on the elastic properties of the macroscopic behavior [35]. The particle mass can be adjusted to verify two conditions; the first is related to the inertial number as follows: in which is the strain rate of the boundary conditions used to simulate the cyclic triaxial testing, and is the volume of the sphere. In an ideal quasi-static condition, 0; Radjai and Dubois [36] suggested that the quasi-static limit is approached for 10 . The second condition is aimed to reduce the computational cost through the upper limit of the time step in order to ensure the stability of the explicit integration scheme through the following relation: where is the particle density linked to the mass, . Therefore, the mass of each particle could be set as high as possible in order to reduce the computational cost of the simulation if the quasi-static condition is met by limiting the inertia number of Equation (15). Thornton and Antony [37] suggested increasing the particle density up to a fictitious 10 12 value, whilst Macaro and Utili [38] adopted a value of 10 9 for the simulation of triaxial tests of seabed sands. On the other hand, the primary parameters are related to the assembly's frictional behavior, such as the interparticle angle of the friction, , and the limiting rolling coefficient, . In particular, the latter has the strongest influence on the peak strength, the residual or critical strength, and the contractive/dilative behavior [20]. These two parameters are calibrated in this paper to define a macro-category of soil whose physical characteristics are consistent with the average properties obtained from real sand with uniform grains. Therefore, the Discrete Element Method is able to capture several complex constitutive behaviors with few parameters; the final behavior will be essentially affected by the initial void ratio and the distribution of the grains of the soil sample, as shown in the following Section 3.

Results
In this section, the results of the cyclic triaxial testing are presented in order to propose a methodology to derive the numerical stiffness and strength degradation relation that can then be used at the mesoscale level to calibrate constitutive models or at a macroscale to calibrate macro-models such as those formulated through the dynamic p-y curve approach [1]. First, monotonic drained triaxial testing was conducted in order to obtain the geotechnical properties of the two investigated samples, and, then, undrained cyclic triaxial testing was performed to derive the degradation curves. The simulation was conducted by the open-source code, YADE [39].

Model Description
In this application, two soil samples featuring the behavior of loose sand and medium-dense sand are modelled through the Discrete Element Method, as shown in Figure  2. The loose packing of the spheres is characterized by a relative density of 33% and a void ratio of 0.72; the medium-dense sample has a relative density of 55% and a void ratio of 0.60. The two soil samples are representative of uniformly graded, well-rounded, coarse sand. Numerical simulations have been conducted on a representative cell with periodic boundary conditions. In order to satisfy periodicity conditions, a periodic space is created by the repetition of a parallelepiped-shaped cell. The use of the periodic boundary allows for the production of homogeneous and isotropic states, eliminating the disturbances in the granular structure induced by the wall effect [12]. The threedimensional assemblies are composed of 20,199 spheres for the loose sand sample and 20,898 for the medium-dense sample, having a mean particle radius of 1 mm (hence, representing coarse sand); the radii, ranging from a minimum of 0.92 mm to a maximum of 1.09 mm, are drawn randomly from a uniform distribution. This approach has been used to avoid the phenomenon of "crystallization," which occurs with the perfectly regular packing of particles with the same radius, yielding an unrealistic, extremely stiff, and strong soil medium. The relative density, , is obtained as follows: , where is the current void ration, and and are the maximum and minimum void ratios, respectively. For uniform spheres with the same radius, 0.90986, and 0.35047 [40]. The current void ratio is expressed as 1 (18) in which the porosity of the sample, , is numerically computed from the subtraction of the total volume of the cell and the total volume of the rigid spheres. The input parameters used in this study are reported in Table 1, defined as described in Section 2.3, and selected in the range of values proposed in [20] to simulate loose and medium-dense sands. The preparation procedure consists of an isotropic compression of a randomly generated packing of floating spheres as follows: (1) initially, the spheres are generated randomly and placed within a cube of a 0.07 m side in order to avoid any contact between them; (2) the packing is subjected to an isotropic compression by moving all the sides (boundary conditions) of the cube uniformly until the target value of the relative density has been reached; (3) the overlaps are eliminated by slightly decreasing the particles' radius by a factor of 0.999 in order to obtain an initial stress-free packing.

Monotonic Triaxial Testing
After the preparation procedure to obtain the target initial relative density, the two specimens are subjected to drained triaxial compression to derive the main geotechnical properties, such as the angle of internal friction and shear modulus. A monotonic triaxial test comprises two phases: (i) isotropic consolidation and (ii) the deviatoric stage. A constant loading strain rate of 0.0001/s is applied to the boundary conditions in order to avoid inertial effects and maintain quasi-static conditions; the inertial number calculated through Equation (15) is equal to 1.65 10 ; hence, smaller than 10 , as aimed. In this study, the timestep used for the analyses is set to 0.01481 s, corresponding to 50% of the limit value defined by Equation (16).
During the isotropic stage, the strain-control conditions are applied (by uniformly moving the artificial sides of the cube) to achieve constant isotropic stress equal to 100 kPa. Subsequently, a deviatoric phase is carried out by moving the upper and bottom periodic boundaries at the constant loading strain rate while controlling the lateral sides in order to maintain the confining stress of 100 kPa. The frictional angle is calculated through the following expression: where q and p are the deviatoric and the mean stress computed at the peak or at the critical state, respectively. Calculated geotechnical parameters have been summarized in Table 2. The results on the strength and stiffness evidence a good agreement with the experimental data for loose and medium coarse-grained sand [41,42]; therefore, the discrete element model is able to capture the global behavior of the soil at the mesoscale. Figure 3 illustrates the stressstrain and volume change behaviors of the two samples in loose and medium-dense states. A maximum value of an axial strain of 20% is applied to reach a critical state. Figure  3a shows that at large strains both loose and medium-dense soils reach the same critical shearing resistance characterized by the residual angle of friction of 27° and the same void ratio, termed the critical void ratio ( 0.75), corresponding to a critical relative density of 29%. This is consistent with the critical state concept [23] in which cohesionless soil tends to move toward the critical void ratio regardless of its initial value of the void ratio. The critical void ratio marks the boundary between a contractive and a dilative response; therefore, it is expected that a very small dilatancy occurs for the loose sand sample ( 0.72) and a much larger dilation for the medium-dense sample ( 0.60). Moreover, it can be noted that the sample at medium-dense density mobilizes a peak shearing resistance of 35°, which is much higher than the residual fraction angle of 27°, whilst an almost constant deviatoric plateau is obtained for the loose soil sample. The evolution of the volumetric strain with the axial strain is shown in Figure 3b; the increment of the volumetric strain indicates an increment of the total volume of the specimen and, hence, a dilative behavior, which is higher in medium-dense sands than in loose sands. For small values of axial strain, a contractive behavior is observed. Figure 4 shows the degradation curve of / versus the shear strain. The tangential shear modulus, , is defined at a strain level of 10 −4 . Typically, the normalized shear modulus, / , decreases with the increase of the shear strain, and it can be observed that the loose sand shows a stronger nonlinearity than the medium-dense sand, as expected. These curves can be used for typical problems of seismic site response analysis where the definition of an equivalent linear model for the soil is required [43].

Cyclic Triaxial Testing
The undrained cyclic triaxial test is carried out by applying a one-way strain-control input, as shown in Figure 5, for the loose and medium-dense sand, respectively; different maximum axial strain amplitudes, , ranging from 0.1% to 0.5% for loose soil and from 0.1% to 5% for medium-dense soil, are applied to identify the effect of the magnitude of the input (related to the geotechnical concept of the cyclic stress ratio) on the degradation as a function of the number of cycles. The same soil specimens of the input parameters in Table 1, tested under monotonic conditions, are now considered for cyclic testing. The first stage of isotropic consolidation is performed under the same drained conditions of the previous test, reaching a confining pressure of 100 kPa. On the other hand, to perform the undrained shearing stage, a constant volume approach is used [26]; instead of simulating the interaction of the solid particles with the fluid occupying the pores of the sample, this approach aims to simulate the main kinematic condition arising from an undrained test, i.e., the shearing deformation occurs, preserving the total volume. Therefore, this condition can be achieved by controlling the prescribed strain of the boundaries in order to maintain the following relation: in which and are the strains of the lateral boundaries of the cell, whilst is the axial strain of the top and bottom boundary controlled during the triaxial testing, according to the one-way input of Figure 5. This method reduces the computational complexity of a solid-fluid couple model, and it is appropriate for simulating the undrained behaviors of granular materials with coarse grains [26], as were the soil samples investigated in this paper. Figure 6a shows the hysteresis curves obtained from the results of the cyclic triaxial testing for the loose sand at several levels of maximum strain. The hysteresis curves show a typical oval-shaped behavior characteristic of loose sand [44]; the area covered by each loop, related to the dissipation of energy, increases with the increase of the level of strain. It can be observed that the hysteresis loops rotate after each cycle, manifesting an important degradation of the soil strength; in real loose sand, this phenomenon is induced by the increment of the excess pore pressure caused by the compaction of the soil grains, which decreases the effective soil stresses. This effect has been properly captured, as shown in Figure 6a, by the numerical simulation conducted under the constant volume condition. To assess the strength degradation, Figure 6b depicts the evolution of the deviatoric stress with the progression of the cyclic test, i.e., the number of cycles. The maximum value of deviatoric stress decreases after each cycle: for a small level of strain above 0.2%, a few cycles are needed to reach a level of liquefaction where the soil loses its strength entirely.
(a) (b) Figure 6. Results of the cyclic triaxial testing for loose sand: (a) deviatoric stress-strain hysteresis curves; (b) evolution of the deviatoric stress with the increase of the number of cycles.
The variation of the deviatoric stress with the number of cycles is due to the re-organization of the granular material and, hence, to the different exchange of contact forces between the particles; the evolution of the normal contact forces for a representative volume element of the loose sand sample is shown in Figure 7. A little variation of the maximum magnitude of the contact forces is observed at a small level of axial strain during the first three cycles; on the other hand, with the increase of the maximum level of strain, the degradation of the maximum value of the contact forces becomes relevant after a few cycles. It can also be observed that the number of interactions between the particles decreases with the increase of the number of cycles, similarly to what occurs in real soil during a liquefaction phenomenon.  Figure 8a illustrates the deviatoric stress-strain cyclic response for medium-dense soil, and Figure 8b shows the evolution of the deviatoric stress with the progress of the cyclic testing (Dr = 55%). It is interesting to note the different shapes of the medium-dense samples compared to the loose ones: they present the typical "S shape" or "banana shape" observed in real experiments [45,46]. The change in the stiffness and strength response with the increase of the strain level corresponds to the tendency to move from the volume contraction to the volume dilation. Remarkably, the samples subjected to low axial strain corresponding to 0.1-0.5% show an important and gradual loss of shear resistance whilst at larger strain levels, the hysteresis loops after an initial partial degradation and become stable, showing a behavior similar to the plastic shakedown of the traditional elastoplastic theory. Moreover, the cyclic test with medium-dense samples reveals a general strain hardening type of behavior: while at low strain levels, the contraction induces a partial degradation, and at a large strain, dilatancy helps to reduce the effect of the cyclic strength degradation.

Discussion
In this section, the results of the cyclic triaxial testing are elaborated to derive the cyclic stiffness and strength degradation curves. Adopting the general fatigue-based approach, these derived curves can be used in models based on the total stress approach where the cyclic degradation can be expressed as a function of the number of loading cycles [47][48][49]. Figure 9a shows, for the loose sand samples, the evolution with the number of cycles of the Strength Degradation Index, , defined as: , where is the maximum deviatoric stress recorded at the first loop and is the maximum deviatoric stress recorded after N strain-controlled loops. It can be observed that the soil strength rapidly decreases with the number of cycles and with the increase of the maximum axial strain, . The authors proposed to consider soil failure when the strength degradation index, , is lower than 0.1. Moreover, it can be observed that the almost linear degrading behavior of the strength before the soil fails. Figure 9b shows the Stiffness Degradation Index, , for the same results on loose sand samples, defined as , (22) where is the initial tangent shear modulus calculated at the first loop, and is the initial tangent shear modulus at the current N loop. Similar to the strength degradation index, the curves show the progressive degradation of the stiffness with the increase of the number of cycles and strain level.  Figure 10a shows the Strength Degradation Index from the results of the mediumdense sand samples. The failure due to the soil liquefaction ( 0.1) occurs at low strain levels when the soil manifests a compressive behavior; this leads to a fast pore pressure buildup and, consequently, a degradation of the soil strength. On the other hand, because of the relevant dilative behavior that occurs at larger strain levels, the degradation of the soil strength is limited, and the samples do not reach a failure condition. This phenomenon affects the tangent stiffness value at each loop; at low strain levels, a softening behavior is observed whilst cyclic hardening is manifested at a large strain with relevant increments, as evidenced in Figure 10b. Finally, the results presented above are elaborated to obtain useful laws for soil damage models. A linear degradation model could be used to capture the relation between the Strength Degradation Index and the number of cycles, as depicted in Figure 11a,b, for loose and medium-dense sand, respectively. Considering the samples in which the failure occurred, the following relation is proposed: where is the current number of cycles, and is an empirical coefficient depending on the maximum level of strain, , which can be expressed through the following power law: in which 2.305 and 1.177 for loose sand, and 0.2241 and 0.59337 for medium-dense sand.

Concluding Remarks
The potential of the discrete element method in simulating the kinematics of an assembly of particles has been exploited in this study to investigate the cyclic behavior of two ideal soil samples representing uniform sand at two different levels of compaction, a loose and a medium-dense state. Few parameters have been used to create the model; the primary ones are the limiting rolling coefficient, which simulates the ability of a particle to roll over the other and, hence, is somehow related to the texture of the soil grains, and the void ratio of the packing; the secondary ones, such as particle density or interparticle stiffness, are fixed according to few computational rules. First, the samples have been tested to characterize the main geotechnical parameters, such as the peak and residual angle of friction, as well as the shear modulus. Second, cyclic triaxial testing was carried out by adopting a constant volume approach to replicate undrained conditions. The degradation of the deviatoric stress and shear modulus has been monitored to obtain a relation with the number of cycles applied to the sample. Briefly: 1. Geotechnical properties such as the peak and residual angle of friction, as well as the shear modulus, are consistent with the values given in the literature from real soils. 2. The different volumetric behavior of loose sand and medium-dense sand has been simulated without resorting to complex constitutive models or the calibration of several parameters. The same input parameters were used for both specimens, evidencing how the distribution of the particles and, hence, the void ratio of the packing is the main aspect affecting the soil behavior of the sands. 3. The cyclic contractive tendency of the loose soil under undrained conditions at small values of strain observed during the monotonic triaxial testing leads to a cyclic degradation of its strength and stiffness; this is compatible with the real phenomenon of the liquefaction or cyclic mobility. 4. The cyclic tendency to dilate in medium-dense sands yields to a stabilization of the soil degradation; this is compatible with the reduction of the excess pore pressure, as occurred in the real medium-dense sands. 5. The coefficients to model a linear cyclic damage model were obtained as a function of the maximum applied axial strain and soil consistency.
Whilst the results presented in this paper are limited to a simplified model of uniform sands, the proposed investigation can be used to simulate a multitude of packing of soils at different states of compaction and dispersion of the grains. Therefore, this approach can be adopted to derive statistical relations between the cyclic soil degradation or hardening and the number of cycles, to help with the calibration of geotechnical models at the mesoor macro-scale.

Conflicts of Interest:
The authors declare no conflict of interest.