Linking the macro-scale response of granular materials during drained cyclic loading to the evolution of micro-structure, contact network and energy components

This study has considered the behaviour of granular materials subjected to drained cyclic loading under constant mean effective stress. Using the discrete element method, cubical, isotropically compressed samples were subjected to 50 loading cycles at different values of mean stress (p′=\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p' =$$\end{document} 100, 200, 300 kPa) and different loading amplitudes (ζ=\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta =$$\end{document} 5%, 10% and 20% of p′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p'$$\end{document}). At low cycle numbers, the deformation mechanism is controlled by contractive volumetric strains, before transitioning to the ratcheting regime, characterised by the persistent accumulation of plastic strains. An energy/work analysis showed that the volumetric work per cycle decreased as hysteresis loops tighten. During ratcheting, most boundary work was dissipated by contact sliding. The mechanical response was controlled by ζ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta$$\end{document}, with little to no influence of p′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p'$$\end{document}. For ζ=5%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta = 5\%$$\end{document}, deformations were confined to the elastic range, with no increase in secant stiffness Gsec\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_{sec}$$\end{document} or shear strength after cyclic loading. For ζ=10%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta = 10\%$$\end{document}, Gsec\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_{sec}$$\end{document} and the shear strength increased after cyclic loading, although no significant expansion of the yield surfaces was observed. The largest loading amplitude (ζ=20%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta = 20\%$$\end{document}) induced yielding at low cycles, leading to significant changes in the fabric, volume and yield surfaces of the samples, and a significant increase of shear strength and Gsec\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$G_{sec}$$\end{document}. At the micro-scale, graph theory was used to quantify the evolution of the contact network. After ∼20\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 20$$\end{document} loading cycles, the network reached a steady-state of constant but persistent topology changes in the material, with most of the topology retained between loading cycles.


Introduction
Cyclic loading must commonly be considered in civil engineering design analysis. While high frequency cyclic loading may be imposed upon soil during earthquakes, lower frequency cyclic loading is experienced adjacent to structures such as wind turbines with monopile foundations and integral bridge abutments, e.g., [5]. The load cycles may vary in amplitude, frequency and total number over the lifetime of the structures. Large amplitude and high frequency loading may cause liquefaction (considered an undrained response), e.g., during earthquakes [15]. Small amplitude and low frequency loading cycles, observed during soil-structure interaction, lead to particle rearrangement and an increase in soil stiffness and pressure [4].
An important factor during drained (low frequency) low amplitude cyclic loading is the accumulation of permanent (plastic) strains or deformation even after a large number of cycles, which is also known as ratcheting ( [2,3,10,13,33]). [33] suggested that, with increasing cyclic load This article is part of topical collection: Energy transport and dissipation in granular systems. amplitudes, the behaviour of granular material will first be purely elastic, where the deformations are fully recoverable, followed by shakedown, where strain accumulates for a finite number of cycles. Moreover, [10] proposed a post-compaction stage where, in the initial cycles, the energy dissipation and the rate at which deformation is accumulated gradually decrease until a low constant value is reached. Then, in the final stage, the accumulation of permanent strain is small but persistent. They call this stage the ratcheting regime.
Ratcheting behaviour has been studied in a number of laboratory studies and has been considered in the development of constitutive models. [18] monitored the change in temperature and horizontal stresses adjacent to integral bridge abutments and was able to capture the strain ratcheting behaviour. The decrease in plastic strains between consecutive cycles was attributed to the densification of the sample. [13] studied small-strain cyclic loading adjacent to wind turbines with monopile foundations. They integrated ratcheting into a plasticity model. This model was modified so that it was able to capture the reduction in the ratcheting rate and tightening of the hysteretic loops related to the dissipation of energy and the increase in secant stiffness.
In this paper cyclic loading at a constant mean effective stress is investigated using the discrete element method (DEM). This study builds upon prior contributions including [25] who compared DEM simulations with cyclic triaxial tests on samples of steel ball bearings to show that DEM can reliably capture the response of granular materials subject to drained cyclic loading. [37] carried out strain-controlled cyclic loading simulations and observed the continuous build up of lateral stresses in loose and dense assemblies. Their simulations were qualitatively comparable to laboratory data produced by [36]. [24] carried out 2D DEM cyclic simulations on gravel-sand mixtures with different fines contents. By studying the contact force anisotropy and the change in strong and weak contacts, they found that the axial stress controlled the strong contacts and the confining pressure had an influence on the weaker contacts.
The dissipation and accumulation of energy can be predicted using continuum mechanics with constitutive models that make an assumption about energy dissipation either implicitly or explicitly. For instance, [31] simulated the energy dissipated due to hysteresis during cyclic loading using IC.G3S, a cyclic, non-linear constitutive model. In DEM, energy tracing is done at the contacts (strain energy and frictional energy) and/or as boundary work which is the energy that is input into the system. In this study, energy tracing is carried out using a modified version of the code granular LAMMPS [28], introduced by [11], to study drained triaxial shearing. [17] demonstrated that the volumetric and distortional components of the boundary work during cyclic loading can be estimated for an undrained (constant volume) system.
Even though DEM is able to capture energy dissipation, few studies have used the method to explore cyclic loading and ratcheting. In this study, cyclic, stress-controlled loads were applied to a cubic sample with periodic boundaries using DEM. The model was capable of replicating the ratcheting behaviour observed in laboratory experiments. In the next section, the numerical model is described, followed by a description of the macro-mechanical response; then, we link the macro-scale response to the micro-scale in terms of the evolution of the coordination number, fabric anisotropy, and the contact force network. The energy components in the system are linked to the macro-and micro-mechanical behaviour of the sample and its contact network.

Numerical model
The initial DEM specimen was created using a cloud of 19,463 non-contacting spheres with a particle size distribution (PSD) similar to Toyoura sand, as shown in Fig.1a. The algorithm used to place the particles in the specimen ensured an homogeneous packing throughout the sample. Periodic boundary conditions were used with a servo-algorithm to control the stress state. The virtual Fig. 1 a Isotropically compressed DEM sample and b particle size distribution for Toyoura sand cubic specimen was isotropically compressed to three different mean effective stresses ( p ′ = 100 kPa, 200 kPa and 300 kPa) with a Poisson's ratio ( ) of 0.2 and a shear modulus (G) of 29.17 GPa. This corresponds to a Young's Modulus of E = 2G(1 + ) = 70 GPa. Similar stiffness parameters were adopted in [22] where a shear modulus of 29.0 GPa was applied and in [32] where E = 70 GPa. The inter-particle friction coefficient ( ) was 0.25. This friction coefficient was selected as it lies in the range between 0.12 and 0.35 for quartz-based sand grains measured by [30]. Additionally, of 0.25 was selected as it is a value that is often used in DEM for silica sand such as in [17]. [14] showed that using a high value of (i.e., ≥ 0.5) can result in a non-physical response. The three initial samples were classified to be in a loose state with an initial void ratio (e) in the range of 0.661 to 0.663. This was confirmed by the observed stress-strain behaviour which is typical of a loose sample. The virtual samples have the same properties as those considered in [14] and the critical state observed in monotonic tests on the samples considered in the current study come very close to the critical state line reported in [14]. This implies that the initial state parameter, 0 , is about 0.032. The contact model for the simulations was the Hertz-Mindlin model based on [12]. Fig.1a is a representation of one of the isotopically compressed samples. The modified version of the DEM code granular LAMMPS [28] introduced by [11] was used.
After isotropic compression, samples were subjected to triaxial cyclic loading at a constant mean effective stress. Each test was carried out with a specific ratio of loading amplitude ( = �amp p � ). This ratio is defined in a similar manner to [34]. However, it is difficult to directly compare the data generated with specific experimental data as we know of no published dataset for constant p ′ , drained, triaxial, two-way cyclic loading. The ratio defines the amplitude of ′ zz applied during the load cycles, while keeping � xx = � yy . For instance, for p ′ =200 kPa and = 10% , the loading amplitude ′amp is equal to ± 20 kPa, and ′ zz varies between 180 and 220 kPa. Four different positions in each cycle were recorded for analysis: (i) the neutral (isotropic) position at the start of a cycle, (ii) the loaded position when � zz = p � + �amp , (iii) the neutral position when the sample is unloaded and (iv) the unloaded position when � zz = p � − �amp . To keep a constant mean effective stress, ′ xx and ′ yy decrease when ′ zz increases and vice-versa, as indicated in Fig. 2.
The amplitudes are indicated in Table 1 and a total of 50 cycles at a constant frequency f = 4 Hz were carried out for each of the 9 simulations. This frequency ensured that the sample remained in a quasi-static state and is identical to the frequency selected in [17] which used the same version of the granular LAMMPS code. The model does not consider rolling resistance, nor does it capture the particle abrasion or breakage, considered by [9] and [7], respectively, that may be observed during cyclic loading.
The inertia number (I), calculated as shown in Eq. 1 according to [27] and [6], remained I < 10 −4 for all the simulations, verifying that quasi-static conditions were preserved during cyclic loading: where ̇ is the strain rate, d 50 is the median particle diameter of the material ( 0.22 mm ), is the particle density ( 2650 kg m −3 ). ̇= zz ∕ t , with zz corresponding to the strain in the z direction during the loading stage, occurring  over t = 0.0625 s , i.e., one quarter of the cycle period of the first cycle.

Ratcheting and hysteresis loops during cyclic loading
In Fig. 3 the accumulation of volumetric strain ( v ) with axial stress ( zz ) can be observed over 50 cycles. The convention here is that an increase in volumetric strain means a contraction of the sample. At a given p ′ the amount of permanent (plastic) volumetric strain increases with increasing . When p ′ is varied while keeping constant, similar levels of v are obtained at the end of cyclic loading. Thus, it can be said that samples with identical exhibit similar behaviour. The increase in v during loading is larger than during unloading, i.e., the sample exhibits a softer response during loading than during unloading. Furthermore, a lower mean effective stress leads to a more rapid increase in v at low cycle numbers ( < 10 ). Despite the fact that the model does not consider particle abrasion or breakage, there is clear evidence of ratcheting behaviour in the simulation data. A decrease in the rate of accumulation of volumetric strain can be observed with an increase in cycle number. The persistent increase of v can be quantified in terms of void ratio. For instance, the initial void ratio (e) for the sample at p � = 300 kPa was 0.663. After 50 cycles the void ratios for the three different loading amplitudes = 5%, 10%, 20% ( ′amp of 15, 30 and 60 kPa) were 0.662, 0.656 and 0.639 respectively. Similar trends were followed by the simulations carried out at p � = 100 and 200 kPa.
The strains were decomposed into their elastic (recoverable) and plastic (permanent) components. The total and plastic strain are indicated in Fig. 4. The plastic strain was obtained by calculating the difference in axial strain between two subsequent cycles. The elastic strain was obtained by calculating the strain that is recovered in each cycle which corresponded to the difference between the total and plastic strain. The elastic strains became virtually  constant after ∼ 20 cycles, in line with the onset of the ratcheting regime observed by [10]. The plastic strains were small in magnitude but persistent, even at high numbers of cycles. As expected, larger strains were observed for the larger loading amplitudes. Still, the magnitude of the strains for = 20% were significantly larger than those observed for = 5% or 10%, which suggests a transition in the underlying deformation mechanism as the loading amplitude increases. Fig. 5c shows the change in loop area. For = 5% and = 10% the change is very small in comparison to the largest amplitude. The loop area plateaus at larger cycle numbers, which reflects the cyclic strain that plateaus at the same time and a reduction in increase in permanent strains. Figure 6 illustrates the hysteresis loops in the deviatoric stress/strain space for one of the simulations ( p � = 200 kPa and = 20% ). The deviatoric stress, q = ±60 kPa in this example, is calculated as defined in [17] to allow the distortional work to become negative: The first cycle is indicated in green and the last cycle is coloured red. As the cyclic loading progresses there is a clear change in the shapes of the hysteresis loops, with a significant decrease in the loop area. These observations are characteristic of ratcheting behaviour [1].
Next, the variation of the secant stiffness G sec during cyclic loading was studied. Previous work carried out by [1] and [19] demonstrated an increase in G sec with cycle number. The secant stiffness G sec for a specific cycle is the slope of the straight line between the point of minimum (negative) deviatoric stress to the point of the maximum (positive) deviatoric stress; see dashed lines in Fig. 6. For of 10% and 20%, G sec increases with increasing cycle number as the hysteretic loops tighten (i.e., become smaller) and progressively become more narrow (Fig. 7). No clear change in G sec is observed for the lowest loading amplitude ( = 5%). This indicates that the induced loading amplitude and subsequent strains were too small to generate significant changes within the sample. The secant stiffness at the end of cyclic loading was directly proportional to p ′ , but inversely related to , with the sample at p � = 100 kPa and = 20% exhibiting the lowest G sec . The difference in secant stiffness in cycle 0 can be observed in Fig. 7a, b, c. The much smaller deviatoric strain for of 5% leads to a much larger initial secant stiffness in comparison to the largest amplitude where the stiffness is lower due to the degradation of stiffness we expect in any granular material subjected to shear deformation.  [16] defined three different yield surfaces called Y 1 , Y 2 and Y 3 . The perfectly linear elastic zone, characterised by very small and fully recoverable strains, is enclosed by the yield surface Y 1 . G sec is constant in this range. The transitional zone between the recoverable and irrecoverable strains is enclosed by Y 2 . In this zone, the material is "non-linear elastic" with hysteretic behaviour, and [16] hypothesised that energy loss can be observed due to small-scale yielding and friction at particle contacts. When Y 2 is engaged the hysteresis loops are prevented from closing and irreversible plastic strains develop [20]. Lastly, Y 3 encloses the zone where more significant plastic strains are developed. It is assumed that in this zone particle sliding is prevalent. However, previous laboratory and FEM analyses by [16] and [13] were not able to assess particle sliding in this zone. In order to identify Y 1 , Y 2 and Y 3 the isotropically compressed DEM samples were sheared at a speed of 0.0001 m s −1 with the same mechanical parameters used during cyclic loading, as shown on Fig. 8. The linear portion of the stress-strain curve was identified for each simulation. This portion marks the elastic zone and is enclosed by the yield surface Y 1 which was at an axial strain ( a ) of 8.5 × 10 −5 , 1.4 × 10 −4 and 1.6 × 10 −4 % for p � = 100, 200 and 300 kPa, respectively. These values are reasonable when compared with the maximum strain range defined by [16], which was between 10 −4 % and 6 × 10 −4 %. Values obtained in this study may be slightly lower due to the use of spherical particles.

Yield surface
During cyclic loading Y 1 was crossed in every cycle for all loading amplitudes. For the two higher loading amplitudes ( = 10% and 20% ), at low cycle numbers ( < 10 ) the hysteresis loops did not close, indicating that Y 2 is engaged. Furthermore, for = 20% , Fig. 9c shows that Y 3 is crossed during the first ∼ 5 cycles which is in agreement with the plastic strains developed during the same cycle range (see Fig. 5). For smaller loading amplitudes, strains developed during cyclic loading are far below the Y 3 yield surface (Fig. 9a, b), remaining mostly within the elastic range of deformation.
Thus, the increase in shear strength with increasing values of loading amplitude can be attributed both to the densification of the material, as well as to the increase in the yield points ( Y 1 , Y 2 , Y 3 ), observed for the largest loading amplitudes. Most of the volumetric strain and stress mobilisation, which result in the increase of shear strength, occurs before the first 10 − 20 cycles. The deformation induced by the smallest loading amplitude was confined mostly to the elastic regime, and for this reason no material strengthening was observed.

Micro-mechanical response
The rose diagrams in Figs. 10 and 11 show the histogram of contact orientations in the material, before and after cyclic loading. Shading of each bin in the histograms corresponds to the average magnitude of normal contact force ( F ave n ) of the contacts orientated within that bin. The contact force in the normal direction for an individual contact is given by: , r a and r b are the particle radii, n is the overlap between particles and n is the unit vector.
After cyclic loading, the orientation of contacts remains homogeneous, typical of isotropic stress conditions, while F ave n decreased for all the simulations tested. The relative decrease in average contact force was largest for samples with = 0.2 , with an average contact force decrease of 6%. The average coordination number (Z) is the average number of contacts that each particle has in the system, and is given by: where N c , the number of contacts, is multiplied by two as each contact is shared between two particles, and N p is the number of particles in the system. The number of particles is constant; thus changes in Z directly reflect changes in the number of contacts in the system. The mechanical coordination number is defined in a similar manner; however it only includes particles that have two or more contacts as it only takes into account the force-transmitting particles. The increase in Z observed in Fig. 12a for the two larger amplitudes suggests that more contacts are created in the system as a result of cyclic loading, with larger increases for higher p ′ and . A slight initial decrease can be observed for the largest which could be related to large plastic strains at initial cycles, leading to particle relocation. This initial decrease can also be related to Fig. 13 where an initial peak in fabric anisotropy can be observed at cycle 1 for the largest amplitude. The increased number of contacts and lower F ave n (observed in Fig. 11) suggests that the stress in the contact network is less centralised, and better distributed among an increased number of contacts between particles. This change in the distribution of the contact forces has been linked to the densification of the material, and to the increase in shear strength observed in Fig. 8. In Fig. 12b the mechanical coordination number follows a similar trend to the coordination number. The changes are slightly larger. These changes in mechanical coordination number are an indication of the change in contact network; however, this variation still remains quite small. [38] clearly show how the fabric anisotropy (A) can influence the mechanical behaviour of granular material. Here, fabric anisotropy, A is calculated as the difference between the maximum and minimum eigenvalues of the second-order fabric tensor Φ pq : where n are unit vectors describing the contact normal orientation for contact j. The value of A approaches zero for isotropic samples. Figure 13 shows the value of A at the neutral position for each cycle. Samples with of 5% and 10% did not exhibit much change in A and remained relatively isotropic. Following an initial sharp increase in A, the sample with the largest loading amplitude became more isotropic after 10 cycles. Recalling that the secant stiffness plateaus after about 10 cycles, it is clear that the largest changes within the sample occurred over those initial 10 cycles.
Particle sliding occurs at contacts when |F t | > |F n | , i.e., when the Coulomb slip criterion is activated. The tangential force F t is calculated incrementally as: where t gives the increment of shear displacement between timestep − 1 and and k t is the contact shear tangential stiffness.  Figure 14 indicates that the percentage of sliding contacts increases initially for of 10% and 20%. This percentage remained relatively constant and below 10% after cycle 20 for all simulations. For the lowest amplitude, the percentage of sliding contacts decreased from the start, reaching values < 1% after ∼ 12 cycles. The small number of sliding contacts is consistent with the linear elastic regime identified in Sect. 3.

Energy dissipation during cyclic loading
The energy input during cyclic loading is given by the boundary work W. The boundary work W at a given timestep is a function of the normal stresses and strains and is calculated incrementally in the DEM simulations as [17]: where the incremental normal strains x , y and z are determined from the changes in position of the periodic boundaries in the model. The boundary work W is decomposed into its two components: the volumetric and distortional work W = W v + W d as proposed by [35]: At a constant mean effective stress ( p ′ ), the change in volumetric strain ( v ) and the change in volumetric work are linearly related according to Eq. 8. The change in distortional work ( W d ) computed from Eq. 9 can also be calculated as the product of the area enclosed by each loading cycle in the qq space ( A loop ) and the sample volume. The deviatoric stress is calculated as described in [17] (by Eq. 2) to allow the distortional work to become negative. Figure 15b shows A loop for an example loading cycle. Figure 15a compares W d calculated from Eq. 9 (solid line) with W d calculated from the area enclosed by each loading cycle (line with markers). A small disparity between the two methods is observed for the first 5 cycles, attributed to the fact that initial cycles  do not close perfectly, leading to minor inaccuracies in the calculation of A loop . In Fig. 16a-c the volumetric and distortional work per cycle decrease with increasing cycle number. The rate of decrease is fastest for the initial 10 cycles. After 20 cycles W v and W d per cycle become relatively constant. This observation agrees with the evolution of elastic and plastic strains (shown in Fig. 5) and with the tightening of the hysteresis loops in Fig. 6. At the largest amplitude ( = 20% ) the volumetric component ( W v ) of the boundary work is initially bigger than the distortional component ( W d ). However, after ∼ 10 cycles, the W d per cycle becomes larger than W v , suggesting that the rate of volume change in the sample decreases at a high number of cycles, and the distortion of the material (captured by W d ) becomes the prevalent deformation mechanism in the system.
The boundary work W, which represents the input of energy into the system, can be either stored inside the sample as strain energy ( E s ) or be dissipated by friction at the sliding contacts E f . Kinetic energy of the particles is negligible. The evolution of these micro-mechanical components of energy, i.e., the recoverable energy stored at the contacts, namely strain energy ( E s ), and the energy dissipated at the sliding contacts, also called frictional energy ( E f ), is studied. The frictional energy is a cumulative term calculated as [11]: 16 W v and W d per cycle at a = 5% b = 10% c = 20% and d Boundary work at p ′ = 300 kPa F t_0 and F t represent the shear forces before and after the Coulomb slip criterion has been applied, respectively. The strain energy ( E s ) is composed of normal ( E sn ) and tangential ( E st ) components. At a given timestep, E sn and E st are given by: Similar to the frictional energy, E st is calculated incrementally. When a contact j between particles is lost, the residual value of E st j is added to E f to preserve the energy balance in the system. It was verified that the error in the energy balance remained well below 1% for all simulations, and therefore, can be regarded as negligible. This minor discrepancy was linked to computational round-off errors. Figure 17a, b shows the evolution of E sn and E st (respectively), while Fig. 17c shows the cumulative frictional energy E f and boundary work W for simulations at p ′ = 300 kPa. Results show that E sn accounts for most of the energy stored at the contacts, while the dissipated frictional energy in Fig. 17c is closely linked to the percentage of sliding contacts (see Fig. 14). Figure 17 shows that most of the energy input into the system (W) is dissipated as frictional energy E f , while the rest is stored as strain energy. Most of this strain energy is stored as normal strain energy ( E sn ). E st decreases with cycle number: consistent with the reduction in sliding contacts.
In Fig. 18 the variation in total stain energy (E s ) with axial strain ( a ) is presented for p ′ of 200 kPa. The same trends were observed for mean effective stresses of 100 and 200 kPa and are omitted to be concise. The total strain energy E s for = 5% increases monotonically during cyclic loading. For = 10% (Fig. 18b), E s increases during the first ∼ 10 cycles, before reaching a plateau, followed by a small decrease at high cycle numbers. The accumulation of axial strain at constant E s is an indication of ratcheting, during which the distortion ( W d ) controls the deformation of the material, and most of the work input is dissipated in the form of E f . Lastly, for = 20% (Fig.18c), both E s and the range of axial strain decrease sharply with cycle number. These results are consistent with Fig. 9, where the decrease in strain energy is associated with initial material yielding, while the reduction in the range of a occurs once the loading cycles fall within the elastic range of the material.

Contact network
The contact network can be described as a directed graph G, with particles and contacts represented as nodes and edges respectively. This analysis presented here focuses on the contact network in the z−direction, following the orientation of ′ zz . The top/bottom periodic boundaries (normal to the z− axis) are represented as virtual nodes in the graph, and particles mapped across the top/bottom boundaries are connected to the corresponding virtual node by virtual edges. The force percolating through the sample in the z−direction, referred as the percolation force f p herein, is the total force transmitted across the sample in the z−direction. Similarly, the ratio between f p and the cross-sectional area of the sample (in the XY plane) is defined as the percolating stress of the network.
The interconnected subset of contacts that transmits f p is defined as the percolation network G perc . G perc is found using the maximum flow algorithm [8], with the top/bottom virtual nodes taken as sink/source, and setting the capacity of the edges to be the z−component of the force at each contact. Thus edges of G that belong to the maximum flow graph (i.e., transmit a non-zero value of f p ), and the value of the maximum flow correspond to G perc and f p respectively. A (a) (b) Fig. 19 Characteristics of the percolation network G perc . a Percolating stress, b Percentage of contacts, and c percentage of particles in G perc detailed description of the implementation of the algorithm is available in [26].
The percolation stress fluctuates in agreement with the imposed variation in zz , as shown in Fig. 19a. Figure 19b, c show that f p is transmitted by 60 − 66% of the contacts in the material, involving 71 − 76% of the particles. The remaining particles and contacts do not contribute to the stress transmission in the z−direction, but may provide lateral support to G perc and contribute to the stress transmission in the x− and y−directions.
The numbers of particles and contacts in G perc are nearly independent of p ′ and , and remain effectively constant during cyclic loading. However, even with a constant number of contacts, there may be an 'exchange' of contacts, with transient contacts leaving/joining the G perc network during cyclic loading. The number of transient contacts is estimated by comparing the contacts in G perc between consecutive loading steps, and the number of new and lost contacts recorded as a percentage of the number of contacts in G perc . New contacts may or may not exist in G before joining G perc , and lost contacts may or may not disappear from G after leaving G perc . Figure 20 shows the variation in the number of transient (new/lost) contacts for the different loading cases. The average number and range of the fluctuation of transient contacts in G perc increases with and is virtually independent of the mean stress ( p ′ ). For the largest loading amplitude  -c), and different loading amplitudes ( ). The initial distribution is common for different values of ( = 20% ), the percentage of transient contacts decreases with cycle number before reaching a stationary behaviour after ∼ 20 cycles, suggesting a steady state is achieved. Conversely, for lower loading amplitudes ( = 5% or 10% ) the percentage (and range of fluctuation) of transient contacts remains constant through cyclic loading. The constant number of contacts in G perc with fluctuating percolation stress (Fig. 19a) suggests differences in the distribution of f p in G perc during cyclic loading. Figure 21 shows the distribution of the percolation force f p at the contacts in G perc before (initial) and after cyclic loading. The distribution of f p is centralised, with a high density of contacts carrying values of f p close to the mean and a small number of contacts carrying low and high values of f p . After cyclic loading, the distribution of f p becomes more uniform, increasing the amount of weak and strong contacts. Changes in the distribution of f p are larger with increasing values of . Figure 22 shows the distribution of f p for the new contacts at different loading stages. The distributions of f p for new and lost contacts are similar, thus lost contacts are not shown. The four different loading steps are: neutral to loaded, loaded to neutral, neutral to unloaded and unloaded to neutral. Neutral state corresponds to the isotropic state, whereas the loaded (+) and unloaded (-) stages correspond to � zz = (1 ± ) ⋅ p � . The mean f p of the transient contacts in G perc increases with . Still, transient contacts predominantly carry small percolation forces. Small changes between loading steps suggest that the mean f p of transient contacts is larger during loading and unloading (leaving isotropic stress conditions) compared to loading steps that return to the isotropic state. Interestingly, the difference between loading steps diminishes with cycle number, once again suggesting a steady state of the fluctuations of the network.
Next, the overall changes in the network were quantified using the dimensionless network change index (NCI), defined according to Eq. 14.
where t is a vector containing the values of f p at every contact in G perc at time step t, t is of constant length and sparse with non-zero entries corresponding to the subset of contacts that belong to G perc at time t.  Fig. 23a, NCI is calculated using the starting points of consecutive cycles. Results are unaffected by p ′ , but significantly influenced by , with higher loading amplitudes inducing larger changes in the network. NCI decreases from an initial value before stabilising after ∼ 20 cycles. These results are consistent with the observed changes in the distribution of f p (Fig. 21) and the percentage of transient contacts (Fig. 20).
For solid lines in Fig. 23b, NCI is calculated in the same way as in Fig. 23a. In contrast, values of NCI in dashed lines in Fig. 23b are calculated between consecutive loading steps (4 loading steps per cycle). Regardless of p ′ and , both the average value and the range of the fluctuations in NCI are significantly larger when calculated between loading steps when compared against NCI calculated at the starting point between cycles. Such results suggest that after an entire loading cycle the networks recover their topology to some extent. In addition, after ∼ 20 cycles, the average value of NCI stabilises, suggesting a recurring but steady level of change in the contact network. Figure 24a shows the distribution of the percentage of transient (new and lost) contacts for each loading stage. Figure 24b shows the evolution of NCI for each loading stage during cyclic loading. At the start of cyclic loading, deviating from the isotropic (neutral) state induces larger changes to the network, both in terms of percentage of transient contacts and NCI. The smallest change to the network occurs when transitioning from the loaded to the neutral state. At larger numbers of cycles ( > 20 ), the differences between loading stages reduces, consistent with the results from Fig. 22. The largest change in G perc , which has a fixed orientation along the z-axis, occurs when the zz becomes the minor principal stress orientation (neutral to unloaded stage). Conversely, at the end of cyclic loading, the stages between the loaded and neutral states induce the least change in the contact network.
The difference between new and lost contacts shows that an increase in zz (unloaded to neutral and neutral to loaded) slightly reduces the number of contacts in G perc while the number of contacts increases when zz decreases. Although counter-intuitive, the results are congruent with Fig. 19c, where the percentage of contacts in G perc decreases with increasing values of . We hypothesise that when the major principal stress of the sample aligns with the orientation of G perc , i.e., at the loaded stage, the percolation network aligns with the principal stress orientation in the material, while the remaining contacts act as lateral support to the principal stress network.

Conclusions
This contribution has considered a series of DEM simulations of triaxial element tests that applied stress-controlled, cyclic loading at constant mean effective stress to a sample of spheres. In the parametric study the mean effective stresses p ′ and the cyclic loading amplitudes were varied. The two main aims of this study were to investigate the macro-and micro-mechanical responses and to further understand energy dissipation during cyclic loading.
The DEM model captured key features of cyclic loading that have been documented in prior experimental studies such as the accumulation of plastic strains, i.e., ratcheting, the rate of which decreases with an increase in cycle number. The hysteresis loops tightened progressively during cyclic loading as their area decreased. This decrease in loop area is related to a reduction in the boundary work per cycle. Most of the boundary work, i.e., the energy input into the system, is dissipated as frictional energy.
The remaining energy is stored at the contacts in the form of strain energy. The use of spheres will lead to a lower deviatoric stress and ' as demonstrated by [21] and [29]. Hence this will also lead to less dissipated energy and a softer response. Other differences might be related to the fact that at the start of cyclic loading all DEM samples are isotropic samples which might not be the case in the lab due to initial fabric. However, the overall trends remain the same. [14] found that particle shape did not affect the shape of the relationship between ' and . The DEM simulations demonstrate that it is useful to interpret the results within the framework of yield surfaces proposed by [16]. For large loading amplitudes ( ≥ 10% in this study), cyclic loading expanded the yield surfaces Y 1 , Y 2 and Y 3 , increasing the elastic range and shear strength of the material. The expansion of the yield surfaces is also linked to the observed reduction in plastic and elastic strains with cycle number.
For samples with the lowest loading amplitude ( = 5%), the coordination number Z and secant stiffness G sec of the material remained virtually constant over the 50 loading cycles, and the fabric anisotropy remained constant and isotropic. Conversely, samples at the highest loading amplitude ( = 20%) exhibited an increase in Z and G sec and a decrease in fabric anisotropy with increasing cycle number.
For = 20% at low cycle numbers ( < 10 ), the volumetric work is the prevalent component of the boundary work. At higher cycle numbers, the deformation is controlled by the distortional work, as elastic strains subside and plastic strains accumulate, i.e., the material starts to exhibit ratcheting behaviour. During ratcheting, most of the energy input is dissipated in the form of frictional energy at the sliding contacts. The loading amplitude controls the evolution of the strain energy ( E s ) in the system. For the two lowest amplitudes ( = 5% and 10% ), E s increased before reaching a plateau at high cycle numbers. For the largest loading amplitude ( = 20% ), E s reached a peak value as the shear strength was mobilised during the first loading cycle, after which it decreased to a plateau.
The evolution of the contact network during cyclic loading shows that the loading amplitude controls the changes in the network with little to no influence of p ′ . After a number of cycles (about 20 in this study) the contact network reached a 'steady state', after which network changes were persistent but constant. At a high number of cycles, the distribution of the percolation force at the contacts, i.e., forces percolating between the top/bottom boundaries of the sample, became more uniform, exhibiting the characteristic behaviour of resilient networks, which could accommodate fluctuations in demand [23].
The changes in the topology of the contact network were determined by the cyclic loading amplitude ; still, regardless of the value of , after 10-20 cycles, the network topology also reached a 'steady state', with persistent but relatively small changes to the arrangements of the microstructure. Interestingly, the contact networks partially restore their topology after a full loading cycle.
The present study unveils the relationship between the macro-scale response of granular media subjected to drained cyclic loading, and the energy and micro-scale mechanisms that ensue in the material. Results show that the amplitude of the cyclic loading determines the response of the material, controlling the changes to its micro-structure, and the subsequent changes in stiffness and shear strength. Regardless of the loading amplitude and at a high number of cycles, drained cyclic loading induces a ratcheting regime, characterised by the prevalent accumulation of plastic strain, distortion of the material, and most boundary work dissipated at the sliding contacts.
Future work may include the development of constitutive models that capture the characteristics and transition into the ratcheting regime, including the energy considerations found at the micro-scale level. Moreover, the stochastic nature of the loading amplitudes encountered in soil-structure applications may motivate the design of resilient granular materials that can be subjected to large number of cycles with adequate levels of stiffness, strength and permanent deformation.
Acknowledgements This work used the HPC facilities at Imperial College London. The first author's project is part of the CDT in Civil and Environmental Engineering at Imperial College London.
Funding This material is based upon work supported by the CDT in Civil and Environmental Engineering at Imperial College. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the CDT in Civil and Environmental Engineering at Imperial College. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising.

Declarations
Compliance with Ethical Standards The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.