3D DEM Simulations of Cyclic Loading-Induced
Densification and Critical State Convergence in
Granular soils

Tao Wang, Antoine Wautier, Chao-Sheng Tang, François Nicot

► To cite this version:

Tao Wang, Antoine Wautier, Chao-Sheng Tang, François Nicot. 3D DEM Simulations of Cyclic Loading-Induced Densification and Critical State Convergence in Granular soils. Computers and Geotechnics, 2024, 173, pp.106559. 10.1016/j.compgeo.2024.106559. hal-04909107

HAL Id: hal-04909107

https://hal.inrae.fr/hal-04909107v1

Submitted on 23 Jan 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1

# 3D DEM Simulations of Cyclic Loading-Induced Densification and Critical State

2

## Convergence in Granular soils

3

Tao WANG<sup>a,*</sup>, Antoine WAUTIER<sup>b</sup>, Chao-Sheng TANG<sup>a,*</sup>, François NICOT<sup>c</sup>

4

<sup>a</sup>School of Earth Sciences and Engineering, Nanjing University, 163 Xianlin Road, Nanjing 210023, China.

5

<sup>b</sup>INRAE, Aix-Marseille University, UR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence, Cedex 5,

6

France.

7

<sup>c</sup>Université Savoie Mont-Blanc, ISTerre, 73376 Le Bourget-du-Lac, France.

8

**Email address:**

9

tao_wang@nju.edu.cn (Tao WANG, Corresponding author)

10

antoine.wautier@inrae.fr (Antoine WAUTIER)

11

tangchaosheng@nju.edu.cn (Chao-Sheng TANG, Co-corresponding author)

12

francois.nicot@univ-smb.fr (François NICOT)

13

14

**Abstract:** Granular soils exhibit very complex responses when subjected to cyclic loading. Understanding the cyclic behavior of such materials is not only crucial for engineering applications but also the bottleneck of most of constitutive models. This study employs 3D Discrete Element Method (DEM) simulations to explore the accumulative plastic deformation and the internal fabric evolution within granular soils during cyclic loading. Two novel observations are identified: 1) A distinct and unique linear relationship between post-cyclic loading void ratio *e* and log (*p**/*p*₀) is found independent of the amplitude of cyclic load and the initial stress state prior to cyclic loading, where *p*\* is the mean pressure incorporating cyclic loading stress and *p*₀ is the mean pressure prior to cyclic loading; 2) When resuming drained triaxial loadings after cyclic loadings, we observe that both microstructural and macroscopic variables converge to the same values they would have reached for pure monotonic drained triaxial loadings. This intriguing behavior underscores and extends to more general loading paths the influential and attractive power of the critical state.

15

16

17

18

19

20

21

22

23

24

25

**Keywords:** Cyclic loading; Discrete Element Method; Accumulative plastic deformation; Anisotropy; Coordination number; Granular material

26

27

28
29

# 1. Introduction

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

Granular soils, including sands and gravels, are subjected to complex cyclic loading in various scenarios, such as waves, wind, earthquakes, construction work or moving traffic (Ng and Dobry, 1994; Donna and Laloui, 2015; Xu and Guo, 2021). Consequently, the study of the cyclic behavior of granular soils has gained considerable attention within the geotechnical community over the past few decades (Lazcano et al., 2020; Song et al., 2023).

Cyclic loads induce the accumulation of plastic deformation in granular soils. Numerous scholars have conducted extensive research on the accumulative deformation of granular soils under cyclic loads. They have explored the effects of various internal influencing factors, e.g., soil type (Li et al., 2021), grain size distribution (Cai et al., 2018), initial density (Xiao et al., 2019) as well as external influencing factors, including cyclic frequency (Thakur et al., 2013), amplitude of cyclic load (Lackenby et al., 2007), principal stress rotation (Guo et al., 2022), and cycle numbers (Indraratna et al., 2012). Among these factors, the amplitude of the cyclic load is highlighted as a critical factor influencing accumulative plastic strain development. The accumulative plastic strain is found to gradually stabilize if the cyclic load is below a critical value (Hyde et al., 1993; Yang et al., 2012). However, when the cyclic load exceeds a certain critical value, it results in a detrimental response of deformation characteristics which makes the specimen collapse (Lei et al., 2016). Furthermore, the soil's initial stress state prior to cyclic loading significantly impacts its dynamic performance. Studies indicate that specimens with higher initial stress ratios or lower confining stresses exhibit greater accumulative axial strain under identical cyclic loading conditions (Peng et al., 2019; Cui et al., 2023). Nevertheless, a comprehensive relationship between cyclic load amplitude, initial stress state, and accumulative deformation of granular soils is still missing. The objective of this study is to develop such a relationship, which holds significant implications for accurately predicting accumulative plastic deformation in granular materials subjected to cyclic loading scenarios.

The mechanisms governing accumulative deformation under cyclic loads in granular materials result from the collective rearrangement of particles through sliding, contact opening, and contact creation in response to cyclic loading (Dean 2005; Kuhn and Chang, 2006; Wautier et al., 2019; Rahman et al., 2021; Mei et al., 2023). In this regard, the discrete element method (DEM), renowned for simulating the movement and interaction of particles within granular assemblies, emerges as an appropriate tool for investigating the micro-mechanical mechanisms driving plastic deformation. Through DEM simulations, Gu et al. (2020) found that the cyclic behaviors of granular materials with the same micro state parameter $\Psi_{MCNO}$, defined as the difference between the initial and critical state mechanical
57

58

59

60

61

62

coordination number, are close to each other, indicating that ΨMCN0 is a plausible state variable for characterizing the
behavior of granular materials. Kolapalli et al. (2023a) performed a large number of DEM cyclic triaxial tests and
evaluated various factors influencing the magnitude and rate of excess pore water pressure generation under cyclic
loading conditions. Additionally, Wang et al. (2021) conducted DEM simulations to investigate the fabric evolution,
plasticity, and dilatancy of sand under cyclic loading and proposed an anisotropic plasticity model based on their find-
ings.

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

.In monotonic loading tests, such as drained triaxial loading (axial compression with constant lateral confining pressure), a granular assembly reaches a critical state—a stationary state where stress and volume tend to remain constant under continuous shear strain. The critical state concept defines a unique linear critical state line (CSL) in e-log (p') space, where e is the void ratio and p' is the effective mean pressure. A critical state soil mechanics framework has proven to be powerful for capturing the monotonic behavior of soils, with some researchers extending it to cyclic loading. These studies use a state parameter (Ψ) for quantifying liquefaction resistance for soil. This state parameter is defined as the difference between the current void ratio and the corresponding void ratio on the CSL at a particular mean effective normal stress (Huang and Chuang, 2011; Zhao and Guo, 2013; Rahman and Sitharam, 2020; Kolapalli et al., 2023b). While the attractor property of the critical state has been studied in monotonic loading (Deng et al., 2021), its application to non-monotonous loading paths, such as cyclic loading, remains unexplored. During stress-controlled cyclic loading, the stress state point (e, p') in e-log (p') space shifts due to accumulative volumetric strain. A key question arises: will a post-cyclic triaxial test drive the material to the same critical state observed in pure monotonic loading? Answering this question is crucial for predicting post-cyclic behavior and establishing constitutive models for granular soils. Investigating the attractor property of the critical state in relation with the significant changes in microstructure induced by cyclic loading is the novel focus of this study.

78

79

80

81

82

83

84

This study intends to explore the accumulative deformation characteristics of granular media under cyclic
loading and its post-cyclic loading behavior by using DEM. The organization of this study is as follows: In Section 2,
drained triaxial tests are first executed, and specific stress states are chosen as initial stress states for subsequent cyclic
loading tests. Section 3 focuses on the analysis of accumulative volumetric strain induced by cyclic loading. In
Section 4, mixed cyclic loading and triaxial loading tests are performed to investigate how a specimen converges
towards its CSL after a cyclic loading test. Finally, concluding remarks are presented in Section 5.
85

# 2. DEM simulation of cyclic loading test

86

## 2.1 DEM model

87

88

89

90

In this study, we utilized the open-source software YADE (Šmilauer et al., 2015) for conducting numerical simulations. The interaction between two grains was modeled using the classical elasto-frictional contact law proposed by Cundall and Strack (1979). The calculation of normal and tangential contact forces ($F_n$ and $F_t$) is outlined as follows:

$$ \left\{ \begin{array}{l} F_n = k_n \delta_n \\ dF_t = k_t d\delta_t, \quad F_t \le F_n \tan\phi \end{array} \right. \qquad (1) $$

91

92

93

where $k_n$ and $k_t$ represent the normal and tangential stiffness, respectively. $\delta_n$ and $\delta_t$ correspond to the relative displacements in the normal and tangential directions, respectively. Additionally, $\phi$ denotes the friction angle that controls the sliding between grains, limiting therefore the tangential contact force through the Coulomb criterion.

94

95

In Equation (1), the normal stiffness $k_n$ is contingent on the size of the two contacting grains, being proportional to a material modulus $E$ and the harmonic average of the radii of the two grains, $R_p$ and $R_q$.

$$ \left\{ \begin{array}{l} k_n = E \frac{2R_p R_q}{R_p + R_q} \\ k_t = r k_n \end{array} \right. \qquad (2) $$

96

97

98

99

100

The parameters employed in this simulation are detailed in Table 1. The parameters selected for this study are in accordance with recommendations given in the Yade software manual (Šmilauer et al., 2015), and are consistent with those used by Xu et al. (2024) and Shi et al. (2024) as well. They were chosen based on a balance between realism and computational efficiency. We did not calibrate the DEM parameters to represent a specific granular material, since our primary objective is to put forward the generic physics of granular materials.

101

<table><caption>**Table 1** Parameters used in DEM simulations</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Density</td><td>3000 kg/m³</td></tr><tr><td>Material modulus (E)</td><td>300 MPa</td></tr><tr><td>Stiffness ratio r = k<sub>t</sub>/k<sub>n</sub></td><td>0.5</td></tr><tr><td>Inter-grain friction angle</td><td>35°</td></tr><tr><td>Grain-wall friction angle</td><td>0°</td></tr></tbody></table>
## 2.2 Sample preparation and cyclic loading test

The grain size distribution for all numerical specimens is illustrated in Fig. 1. The ratio of maximum diameter to minimum diameter is five. To prepare the numerical samples, we initially generate a cloud of 12,000 non-overlapping spheres within a box surrounded by frictionless walls. Subsequently, a consolidation process is applied to isotropically compress the specimen. During the consolidation process, the contact friction is adjusted to control the final density of the samples, wherein smaller friction angles result in denser samples. In this study, we consider loose specimens with a friction angle of 35°, consistent with other references (Lobo-Guerrero and Vallejo, 2006; Wang et al., 2021). Void ratios of numerical specimens after consolidation under different confining pressures are given in Table 2. Then, a vertical compression ($\epsilon_1$) is applied at a strain rate of 0.01/s while maintaining lateral stresses ($\sigma_2=\sigma_3$) constant. Note that the inertial number is below $10^{-4}$ during the whole process of triaxial loading, which ensures that the loading can be regarded as quasi-static (Anandarajah, 2008; Martin et al., 2020). Soil mechanics conventions are adopted where compressions and contractions are counted positive.

<table><caption>**Table 2** Void ratios of numerical specimens after consolidation under different confining pressures</caption><thead><tr><td>Confining pressure</td><td>Void ratio</td></tr></thead><tbody><tr><td>100 kPa</td><td>0.751</td></tr><tr><td>200 kPa</td><td>0.745</td></tr><tr><td>400 kPa</td><td>0.732</td></tr><tr><td>800 kPa</td><td>0.710</td></tr></tbody></table>

115

116
<figure><img src="image_2.png" /><figcaption>Fig. 1 Grain size distribution of 3D numerical specimens</figcaption></figure>

Figure 2 illustrates the progression of deviatoric stress ($q=\sigma_1 - \sigma_3$) and volumetric strain ($\epsilon_v = \epsilon_1 + \epsilon_2 + \epsilon_3$) with axial strain for the numerical specimens subjected to various confining pressures during triaxial loading. It is clear that both $q$ and $\epsilon_v$ consistently rise with axial strain across all curves and demonstrate greater strength and higher contractancy under increased confining pressures. To investigate the response of a specific specimen under cyclic loading at distinct mechanical stress states, different specimens were saved during the triaxial loading, each characterized by a different stress ratio ($\eta=p'/q$). The respective stress states are denoted by solid red circles in Fig. 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2 Deviatoric stress $q$ (a) and volumetric strain $\epsilon_v$ (b) versus axial strain for the numerical specimens in triaxial loading test under various confining pressures. Solid red circles indicate stress states saved for subsequent cyclic loading test.</figcaption></figure>
For the stress states selected during triaxial loading, a cyclic loading test is subsequently conducted in direction **e**₁ as defined in Fig. 1 while maintaining constant confining pressures in direction **e**₂ and **e**₃, as shown in Fig. 3(a). During the cyclic loading process, a constant downward (and upward) velocity is initially applied to the upper (and bottom) wall to initiate the loading stage and impose an incremental stress Δσ₁. Once the stress on the upper wall reaches the target value σ₁+Δσ₁, the movement direction of the upper (and bottom) wall is reversed to initiate the unloading stage. When the stress reduces to the initial axial stress σ₁, the upper wall moves downward again. This completes one cyclic loading cycle, during which the lateral confining stress is kept constant. Various combinations of stress states, characterized by σ₃, σ₁ and Δσ₁, are considered in this study, as outlined in Table 3. Throughout both the loading and unloading stages, the velocity of the upper wall is maintained constant and low enough, corresponding to an axial strain rate of 0.01/s, to ensure quasi-static conditions throughout the entire cyclic loading process. Fig. 3(b) give the axial stress–strain curve during the first 5 cycles for illustration under the case of σ₃=100 kPa, η=0.39 and Δσ₁=10 kPa. It shows that accumulative axial strain greatly develops during the first two cycles and gradually slows down in the following cycles.

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3 (a) Movement of the walls to conduct cyclic loading and (b) Axial stress–strain curve during the first 5 cycles for illustration (σ₃=100 kPa, η=0.39, Δσ₁=10 kPa).</figcaption></figure>

144
145
146
147
148
149
150
151
152
153
154
155
156
<table><caption>**Table 3 Programs of cyclic loading test**</caption><thead><tr><th>σ<sub>3</sub> (kPa)</th><th>σ<sub>1</sub> (kPa)</th><th>η</th><th>Δσ<sub>1</sub> (kPa)</th></tr></thead><tbody><tr><td rowspan="2">100</td><td>145</td><td>0.39</td><td>10, 20, 30, 40, 50, 60</td></tr><tr><td>159</td><td>0.49</td><td>5, 10, 20, 30, 40, 50</td></tr><tr><td rowspan="2">200</td><td>302</td><td>0.44</td><td>20, 40, 60, 80, 100</td></tr><tr><td>343</td><td>0.58</td><td>10, 20, 30, 40</td></tr><tr><td rowspan="2">400</td><td>623</td><td>0.47</td><td>20, 50, 100, 150, 200</td></tr><tr><td>693</td><td>0.59</td><td>40, 80, 120, 160</td></tr><tr><td rowspan="2">800</td><td>1191</td><td>0.42</td><td>50, 100, 200, 300, 400</td></tr><tr><td>1386</td><td>0.59</td><td>40, 80, 120, 160</td></tr></tbody></table>

## 3. Macroscopic volumetric response to cyclic loading

### 3.1 Accumulative volumetric strain

Figure 4 illustrates the evolution of accumulative volumetric strain ($\varepsilon_v^{ac}$) with the number of cycles (N) at various loading amplitudes ($\Delta\sigma_1$) for specimens under different initial stress states (i.e., different $\sigma_3$ and $\eta$). All $\varepsilon_v^{ac}$-N curves exhibit a monotonous increasing and eventually reach a steady value as N becomes relatively large. Note that in all simulations, the values of $\sigma_1+\Delta\sigma_1$ do not exceed the peak stress of the specimen which is obtained at the critical state during monotonic triaxial loading, so that no detrimental response of deformation is observed. With the increase in $\Delta\sigma_1$, both the number of cycles required to achieve a stabilized regime and the corresponding value of $\varepsilon_v^{ac}$ increase. For smaller $\Delta\sigma_1$ (e.g., $\sigma_3=800$ kPa, $\eta=0.42$ and $\Delta\sigma_1=50$ kPa), $\varepsilon_v^{ac}$ stabilizes at a low value of 0.001 when the number of cycles is only 100. However, for large $\Delta\sigma_1$, the resulting accumulative volumetric strain is significant. Considering the case of $\sigma_3=100$ kPa, $\eta=0.39$ as an example, $\varepsilon_v^{ac}$ exceeds 0.03 when $N > 2500$ at $\Delta\sigma_1=60$ kPa. This means that a significant change in microstructure occurs during cyclic loading, consistent with the increase in the accumulative volumetric strain.
<figure><img src="image_6.png" /><figcaption>(a) $\sigma_3=100 \text{ kPa}, \eta=0.39$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) $\sigma_3=100 \text{ kPa}, \eta=0.49$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) $\sigma_3=200 \text{ kPa}, \eta=0.44$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) $\sigma_3=200 \text{ kPa}, \eta=0.58$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e) $\sigma_3=400 \text{ kPa}, \eta=0.47$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(f) $\sigma_3=400 \text{ kPa}, \eta=0.59$</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 4** Development of accumulative volumetric strain with number of cycles for specimens at various loading amplitudes and initial stress states.</figcaption></figure>

## 3.2 Analysis of volumetric response in e-p\* plane

The loose specimens investigated in this study undergo volumetric contraction during the cyclic loading process, resulting in a decrease in the void ratio. The void ratios after the cyclic loading test are computed under various loading amplitudes and are presented in Fig. 5 in terms of p\*. Here, p\* denotes the maximum mean pressure reached by the specimen during the cyclic loading by incorporating the cyclic loading amplitude (Δσ₁), expressed as p\*=(2σ₃+σ₁ + Δσ₁)/3, where σ₁ is the axial stress prior to the cyclic loading test. It is worth noting that the post-cyclic loading void ratios demonstrate a linear relationship with p\* in the e-log (p\*) plane, with a high correlation coefficient R². Indeed, this linear relationship holds true irrespective of the initial stress state, encompassing different σ₃ and η. Similarly, Huang et al. (2019) observed that the final steady value of the void ratio (i.e., minimum void ratio) under cyclic loading varies linearly with the corresponding normal stress both in the shear band and outside the shear band of the specimen. They also noted that the slopes of these linear relationships are significantly different. Henceforth, we shall refer to this line in the e-log (p\*) plane as the maximum densification lines for cyclic loading, since the void ratio reach its minimum value under a certain cyclic load.
<figure><img src="image_14.png" /><figcaption>(a) $\sigma_3=100\text{kPa}, \eta=0.39$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) $\sigma_3=100\text{kPa}, \eta=0.49$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(c) $\sigma_3=200\text{kPa}, \eta=0.44$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(d) $\sigma_3=200\text{kPa}, \eta=0.58$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(e) $\sigma_3=400\text{kPa}, \eta=0.47$</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(f) $\sigma_3=400\text{kPa}, \eta=0.59$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(g) $\sigma_3=800\text{kPa}, \eta=0.42$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(h) $\sigma_3=800\text{kPa}, \eta=0.59$</figcaption></figure>

**Fig. 5** Evolution of the void ratios after cyclic loading at various loading amplitudes as a function of $p^*$. $p^*$ is the maximum mean pressure reached by the specimen during the cyclic loading.

Figure 6 presents all the post-cyclic loading void ratio data in the e-log ($p^*$) plane. A notable observation from Fig. 6 is that the post-cyclic loading void ratios of specimens subjected to the same confining pressure $\sigma_3$ lie on a single line in the e-log ($p^*$) plane. This indicates that after cyclic loading, specimens with an initial stress state corresponding to the same confining pressure share the same maximum densification line. Moreover, the slopes of the four fitted lines (i.e., $\sigma_3=100, 200, 400, \text{and } 800\text{ kPa}$) are close, resulting in nearly parallel lines in the e-log ($p^*$) planes.

<figure><img src="image_22.png" /><figcaption>Fig. 6 Void ratios of specimens after cyclic loading versus $p^*$ under various loading amplitudes and initial stress states.</figcaption></figure>
The data points representing the post-cyclic loading void ratios for all specimens are plotted in the e-log ($\frac{p^*}{p_0}$) plane, as shown in Fig. 7. In this figure, $p^*$ is normalized by its initial mean pressure $p_0$, $p_0 = \frac{2\sigma_3+\sigma_1}{3}$. An important finding is that these post-cyclic loading void ratios consistently form a single, unique line in the e-log ($\frac{p^*}{p_0}$) plane, independent of the initial stress state. This implies that for any stress states during a triaxial loading path, if a drained cyclic loading path is imposed at those stress states, they will all fall on the same maximum densification line. With this finding, by knowing the initial stress state before cyclic loading ($p_0$) and cyclic loading amplitude ($\Delta\sigma_1$), we can directly predict the void ratio after a large number of cyclic cycles.

<figure><img src="image_23.png" /><figcaption>Fig. 7 Void ratios of specimen after cyclic loading versus $p^*/p_0$ under various initial stress states. $p_0$ is the mean pressure of the specimen prior to cyclic loading test, $p_0 = \frac{2\sigma_3+\sigma_1}{3}$.</figcaption></figure>

## 4. Mixed cyclic loading and triaxial loading paths

This section aims to investigate how the cyclic loading path influences a specimen's progression towards the critical state. A fundamental question arises: if a specimen undergoes a cyclic loading path at a certain stress state, will it eventually reach the same critical state, just as it would without experiencing cyclic loading? In other words, does critical state depend on past loading history? To address this question, a comprehensive loading path combining cyclic loading path and triaxial path was simulated, as depicted in Fig. 8. The simulation was conducted as follows: after the cyclic loading paths were applied to the specimens as described in Section 2, a secondary triaxial loading path was executed while maintaining the lateral stress at its current level. These combined loading paths enable to evaluate whether the specimen converges to the same critical state when cyclic loading is applied.
<figure><img src="image_24.png" /></figure>

212

213

214

215

**Fig. 8** Complex loading path combining cyclic loading path and triaxial loading path. The solid blue point denotes the state saved from triaxial loading to initiate cyclic loading. On the other hand, the solid red point represents the state when cyclic loading is completed and serves as the starting point for the subsequent triaxial loading.

216

**4.1. p-q-e analysis**

217

218

219

220

221

222

223

224

225

226

227

The evolutions of deviatoric stress $q$ with axial strain during the process of secondary triaxial loading are depicted in Fig. 9. To facilitate comparison, the $q$ curves from the primary triaxial loading are also included, as shown in Fig. 2 (section 2). One notable observation is that after cyclic loading, the $q$ curves deviate from their original stress-strain paths, exhibiting a sharp increase at the beginning of loading and displaying considerably higher stiffness. This trend aligns with the finding that significant microstructure changes occur during cyclic loading, resulting in denser specimens. Furthermore, it is evident that the deviatoric stress increases with the rise of $\Delta\sigma_1$, some specimens even display stress softening under high $\Delta\sigma_1$, characteristic behavior of “dense specimen”. For instance, cases such as $\sigma_3=100$ kPa, $\eta=0.39$ and $\Delta\sigma_1=60$ kpa, as well as $\sigma_3=200$ kPa, $\eta=0.44$ and $\Delta\sigma_1=100$ kpa exemplify this behavior. Similar findings were reported by Cui et al. (2019), who conducted numerical monotonic simple shear tests on these samples after 6000 cycles of symmetric loadings. They observed that loose specimens densified during the cyclic loading, with the void ratio converging toward those of dense samples.

228

229

230

231

An intriguing observation is that all curves converge to the same deviatoric stress, regardless of the imposed $\Delta\sigma_1$ during cyclic loading history. This deviatoric stress corresponds to the critical state stress obtained during standard monotonic triaxial loading path. It shows that critical state acts as a compelling attractor, pulling the specimen back towards the final state it should reach, despite the cyclic loading path induces significant alterations in the current state,
<figure><img src="image_25.png" /><figcaption>**Fig. 9** Evolution of deviatoric stress with axial strain during the process of secondary triaxial loading after cyclic loading. Curves of reference triaxial loading are included for comparison. Red solid points represent the stress state after cyclic loading (i.e., the beginning of secondary triaxial loading).</figcaption></figure>

Figure 10 illustrates the progression of void ratios for specimens during the secondary triaxial loading process. The critical state line obtained by fitting the critical void ratios under different confining pressures from monotonic triaxial loading is depicted as red curves. Hollow markers are plotted to represent the void ratio at the beginning of the secondary triaxial loading (i.e. end of cyclic loading). Remarkably, all post-cyclic loading void ratios evolve towards the critical state line, regardless of their initial state, and ultimately converge to the same critical state point as achieved in the loading path of monotonic triaxial loading. This feature confirms the remarkable attracting power of the critical state.
<figure><img src="image_26.png" /><figcaption>**Fig. 10** Development of void ratios during secondary triaxial loading. Hollow markers are plotted to show the void ratio of the state at the beginning of secondary cyclic loading (i.e., end of cyclic loading).</figcaption></figure>

The findings reported in Fig. 9 and Fig. 10 provide compelling evidence that specimens evolve invariably towards the same critical state *p*-q-e lines, irrespective of the altering effect of cyclic loading which induces significant fabric changes within the specimen.

**4.2 Microscopic evolution**

The coordination number (CN), introduced by Rothenburg and Bathurst (1989), is a micromechanical metric that quantifies the contact density within granular assemblies. It is defined as $CN = 2N_c/N_{total}$, where $N_c$ represents the total number of contacts and $N_{total}$ represents the total number of particles. The evolution of CN during the secondary triaxial loading under various initial states is presented in Fig. 11. Additionally, the evolution of CN within specimens that only undergo triaxial loading is displayed for comparison marked as black curves. During triaxial loading with
$\sigma_3=100$ and 200 kPa, CN remains below 3 throughout the entire loading process. As $\sigma_3$ increases to 400 and 800 kPa, CN increases overall, indicating that higher confining pressure enhances the contact density within the granular material. After cyclic loading, there is an increase in CN. Subsequently, with the development of axial strain, CN reaches the same value as the specimen during the reference triaxial loading. This observation is consistent with the behavior of $q$ and $e$, both of which converging to the critical state value of the monotonic triaxial loading path. It should be noted that for the sake of clarity, only selected cases are presented here. The results for other cases are aligned with the above findings but not shown here.

<figure><img src="image_27.png" /><figcaption>(a) $\sigma_3=100$ kPa, $\eta=0.39$</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(b) $\sigma_3=200$ kPa, $\eta=0.58$</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(c) $\sigma_3=400$ kPa, $\eta=0.47$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(d) $\sigma_3=800$ kPa, $\eta=0.42$</figcaption></figure>

**Fig. 11** Changes of coordination number (CN) with axial strain during the process of secondary triaxial loading. Hollow circle represents the CN at the state saved for cyclic loading. Hollow square and triangle represent the CN at the beginning of secondary cyclic loading.
Anisotropy is an important concept in granular materials which reflects not only the fabric composition in connection to the spatial arrangement of particles, voids and interparticle contacts but also the changes of these microstructures induced by applied loads (Oda, 1982; Hoque and Tatsuoka, 1998). Within a granular assembly, two main types of anisotropy are recognized: geometrical anisotropy and mechanical anisotropy (Cambou et al., 2004; Rothenburg and Bathurst, 1989). Geometrical anisotropy is defined as the local orientation of a contact plane that gives rise to the global anisotropic phenomenon. It can be quantified by a scalar $a_c$, which represents the deviatoric invariants of the deviatoric part of fabric tensor as introduced by Oda (1982). Mechanical anisotropy is mainly caused by external forces and depends on the induced contact forces in relation to contact plane orientations. It can be quantified by a scalar $a_n$, which is associated with distribution of normal contact force<fnref n="1" />. Well-established formulas for $a_c$ and $a_n$ can be found in Guo and Zhao (2013) and are reviewed in Appendix A1.

Figure 12 shows the evolutions of $a_c$ and $a_n$ during secondary triaxial loading. It can be seen that after cyclic loading, $a_c$ experiences significant increase for loose specimens considered in this study, while $a_n$ keeps almost the same value. This highlights that the cyclic loading process changes the geometrical fabric (characterized by the distribution of contact normal) but have little effect on the mechanical fabric (characterized by the distribution of normal contact force) of the specimen. When a specimen is under cyclic loading, plastic deformation largely depends on the rearrangement of grains by sliding, contact opening, and contact creation, all being related to the geometrical fabric. Therefore, $a_c$ show important change after cyclic loading. However, $a_n$ is more related to the mechanical state of the specimens. The fact that the cyclic tests are stress-controlled is consistent with that $a_n$ does not change after cyclic loading. During secondary triaxial loading, both $a_c$ and $a_n$ evolve towards the critical values of the specimen along the reference triaxial loading path.

<footnote n="1">¹ Note that mechanical anisotropy should involve both normal contact force and tangential contact force. The anisotropy of normal contact force has been proven to be much more pronounced than that of tangential contact force. Thus, we only consider anisotropy of normal contact force here.</footnote>
<figure><img src="image_31.png" /><figcaption>(a) $a_n$ ($\sigma_3$=400 kPa $\eta$=0.59)</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(b) $a_c$ ($\sigma_3$=400 kPa $\eta$=0.59)</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(c) $a_n$ ($\sigma_3$=800 kPa $\eta$=0.42)</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>(d) $a_c$ ($\sigma_3$=800 kPa $\eta$=0.42)</figcaption></figure>

**Fig. 12** Evolution of anisotropy ($a_n$ and $a_c$) with axial strain during the process of secondary triaxial loading. Hollow circle represents the anisotropy at the state saved for cyclic loading. Hollow square and triangle represent the anisotropy at the beginning of secondary cyclic loading.

**5. Concluding remarks**

This manuscript intends to shed the light on the effect of initial stress state and cyclic amplitudes on the accumulative plastic deformation and the internal fabric evolution within granular materials during cyclic loading by means of 3D DEM simulations. The major novelties of this contribution are twofold:

* One significant finding of this research is the identification of a distinct linear relationship between post-cyclic loading void ratio ($e$) and the mean pressure of the specimen incorporating cyclic loading pressure ($p^*$) in the $e$-log($p^*$) plane across diverse cyclic loading amplitude. Furthermore, the post-cyclic loading void ratios were observed to
consistently fall on a unique line in the e-log (p*/p₀) plane, regardless of the initial stress state p₀ prior to cyclic loading, showing the existence of a unique maximum densification line of cyclic loading path. With this finding, by knowing the initial stress state before cyclic loading (p₀) and cyclic loading amplitude (Δσ₁), we can directly predict the void ratio after a large number of cyclic cycles.

• The study also delved into the interaction of cyclic loading with subsequent triaxial loading, revealing that specimens, despite significant microstructure alterations are induced by cyclic loading, eventually converge to the same critical state as they would attain along a unique triaxial loading without experiencing cyclic loading. This intriguing behavior underscores the influential and attractive power of the critical state, pulling the specimen back towards its anticipated critical state despite perturbations from cyclic loading.

In summary, this study illuminates the presence of a unique maximum densification line in granular materials under cyclic loading, providing essential insights for the development of constitutive models and accurate prediction of accumulative plastic deformation in granular materials subjected to cyclic loading scenarios. However, it's important to note that this study primarily focuses on stress-controlled drained cyclic loading paths. Future research should expand the investigation to include other cyclic loading modes, such as undrained cyclic loading and constant *q* (or *p*) cyclic loading, to further enrich our understanding of the complex mechanical responses of granular materials in various loading conditions.

309

## 310 Acknowledgement

This work was supported by the Natural Science Foundation of Jiangsu Province (No. BK20230954), the Open Research Fund of Key Laboratory of Construction and Safety of Water Engineering of the Ministry of Water Resources, China Institute of Water Resources and Hydropower Research (No. 202201), the Research Center for Levee Safety Disaster Prevention, Ministry of Water Resources (No. LSDP202303), and the key laboratory of failure mechanism and safety control techniques of earth-rock dam of the ministry of water resources P. R. China (No. YK323004). For the purpose of Open Access, a CC-BY public copyright licence has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission.

## 319 Reference

[1] Anandarajah A.Multi-mechanism anisotropic model for granular materials. International Journal of Plasticity, 2008, 24:804-846
[2] Cai Y, Chen J, Cao Z, et al. Influence of grain gradation on permanent strain of unbound granular materials under low confining pressure and high-cycle loading[J]. International Journal of Geomechanics, 2018, 18(3): 04017156.

[3] Cambou B, Dubujet P, Nouguier-Lehon C. Anisotropy in granular materials at different scale[J]. Mechanics of materials, 2004, 36:1185-1194.

[4] Cui K, Zhang D, Li P, et al. Influence of the initial static stress state on the accumulation behaviour of a coarse-grained soil under long-term cyclic loading[J]. Soil Dynamics and Earthquake Engineering, 2023, 172: 108042.

[5] Cui, L., Bhattacharya, S., Nikitas, G., & Bhat, A. Macro-and micro-mechanics of granular soil in asymmetric cyclic loadings en-
countered by offshore wind turbine foundations. Granular Matter, 2019, 21, 1-20.

[6] Cundall PA, Strack OD. A discrete numerical model for granular assemblies[J]. Geotechnique, 1979, 29:47-65

[7] Dean, E. T. R. Patterns, fabric, anisotropy, and soil elasto-plasticity[J]. International journal of plasticity, 2005, 21, 513-571.

[8] Deng N, Wautier A, Thiery Y, et al. On the attraction power of critical state in granular materials[J]. Journal of the Mechanics and Physics of Solids, 2021, 149: 104300.

[9] Di Donna A, Laloui L. Response of soil subjected to thermal cyclic loading: experimental and constitutive study[J]. Engineering Geology, 2015, 190: 65-76.

[10] Gu, X.Q, Zhang, J., & Huang, X. DEM analysis of monotonic and cyclic behaviors of sand based on critical state soil mechanics framework. Computers and Geotechnics, 2020, 128, 103787.

[11] Guo N, Yang F, Yang Z X, et al. Deformation characteristics of inherently anisotropic granular media under repeated traffic loading: a DEM study[J]. Acta Geotechnica, 2022, 17(8): 3377-3395.

[12] Guo N, Zhao J. The signature of shear-induced anisotropy in granular media[J]. Computers and Geotechnics, 2013, 47:1-15.

[13] Hoque E, Tatsuoka F. Anisotropy in elastic deformation of granular materials. Soils and foundations, 1998, 38:163-179.

[14] Huang, A.B. and Chuang, S.Y. Correlating cyclic strength with fines contents through state parameters. Soils and foundations, 2011, 51(6), 991-1001.

[15] Huang, M., Chen, Y., & Gu, X. Discrete element modeling of soil-structure interface behavior under cyclic loading. Computers and Geotechnics, 2019, 107, 14-24.

[16] Hyde A F L, Yasuhara K, Hirao K. Stability criteria for marine clay under one-way cyclic loading[J]. Journal of geotechnical engineering, 1993, 119(11): 1771-1789.

[17] Indraratna B, Thakur P K, Vinod J S, et al. Semiempirical cyclic densification model for ballast incorporating particle breakage[J]. International Journal of Geomechanics, 2012, 12(3): 260-271.

[18] Kolapalli, R., Rahman, M. M., Karim, M. R., & Nguyen, H. B. K. A DEM investigation on the influence of cyclic and static stress ratios and state variables on the pore water pressure generation in granular materials. Journal of Geotechnical and Geoenvironmental Engineering, 2023, 149(7), 04023043.

[19] Kolapalli, R., Rahman, M. M., Karim, M. R., & Nguyen, H. B. K. The failure modes of granular material in undrained cyclic load-
ing: A critical state approach using DEM. Acta Geotechnica, 2023b, 18(6), 2945-2970.

[20] Kuhn M R, Chang C S. Stability, bifurcation, and softening in discrete systems: a conceptual approach for granular materials[J]. International Journal of Solids and Structures, 2006, 43(20): 6026-6051.

[21] Lackenby J, Indraratna B, McDowell G, et al. Effect of confining pressure on ballast degradation and deformation under cyclic triaxial loading[J]. Géotechnique, 2007, 57(6): 527-536.

[22] Lazcano D R P, Aires R G, Nieto H P. Bearing capacity of shallow foundation under cyclic load on cohesive soil[J]. Computers and Geotechnics, 2020, 123: 103556.

[23] Lei H, Lu H, Wang X, et al. Changes in soil micro-structure for natural soft clay under accelerated creep condition[J]. Marine Georesources & Geotechnology, 2016, 34(4): 365-375.

[24] Li Y, Nie R, Yue Z, et al. Dynamic behaviors of fine-grained subgrade soil under single-stage and multi-stage intermittent cyclic loading: Permanent deformation and its prediction model[J]. Soil Dynamics and Earthquake Engineering, 2021, 142: 106548.
364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

[25] Li X S, Dafalias Y F, Wang Z L. State-dependant dilatancy in critical-state constitutive modelling of sand[J]. Canadian Geotechnical Journal, 1999, 36(4): 599-611.

[26] Lobo-Guerrero S, Vallejo LE. Discrete element method analysis of railtrack ballast degradation during cyclic loading. Granular Matter, 2006, 8:195.

[27] Martin, E.L., Thornton, C. and Utili, S. Micromechanical investigation of liquefaction of granular media by cyclic 3D DEM tests. Géotechnique, 2020, 70(10): 906-915.

[28] Mei J, Ma G, Tang L, et al. Spatial clustering of microscopic dynamics governs the slip avalanche of sheared granular materials[J]. International Journal of Plasticity, 2023, 163: 103570.

[29] Ng T T, Dobry R. Numerical simulations of monotonic and cyclic loading of granular soil[J]. Journal of Geotechnical Engineering, 1994, 120(2): 388-403.

[30] Oda M. Fabric tensor for discontinuous geological materials. Soils and Foundations, 1982, 22:6-108.

[31] Peng K, Zhou J, Zou Q, et al. Deformation characteristics of sandstones during cyclic loading and unloading with varying lower limits of stress under different confining pressures[J]. International Journal of Fatigue, 2019, 127: 82-100.

[32] Rahman, M. M., Nguyen, H. B. K., Fourie, A. B., & Kuhn, M. R. Critical state soil mechanics for cyclic liquefaction and postliquefaction behavior: DEM study. Journal of Geotechnical and Geoenvironmental Engineering, 2021, 147(2), 04020166.

[33] Rahman, M.M. and Sitharam, T.G. Cyclic liquefaction screening of sand with non-plastic fines: Critical state approach. Geoscience Frontiers, 2020, 11(2), 429-438.

[34] Rothenburg, L. & Bathurst, R. J. Analytical study of induced anisotropy in idealized granular materials. Géotechnique, 1989, 39(4): 601-614.

[35] Shi, X., He, Z., Zhao, J. and Liu, J. Determination of the size of representative volume element for gap-graded granular materials. Powder Technology, 2024, 119578.

[36] Šmilauer V et al. Yade Documentation 2nd ed. 2015, The Yade Project doi:10.5281/zenodo.34073.

[37] Song Z, Yang Z, Zhang M, et al. Numerical characterization of sandstone's mechanical responses under multi-level compressive differential cyclic loading (DCL): A 3D particle-based numerical investigation[J]. Computers and Geotechnics, 2023, 154: 105109.

[38] Thakur P K, Vinod J S, Indraratna B. Effect of confining pressure and frequency on the deformation of ballast[J]. Géotechnique, 2013, 63(9): 786-790.

[39] Thornton, C.. Numerical simulations of deviatoric shear deformation of granular media. Géotechnique, 2000, 50 (1): 43-53.

[40] Wang, R., Cao, W., Xue, L., & Zhang, J. M. An anisotropic plasticity model incorporating fabric evolution for monotonic and cyclic behavior of sand. Acta Geotechnica, 2021, 16, 43-65.

[41] Wang T, Wautier A, Liu S, et al. How fines content affects granular plasticity of under-filled binary mixtures[J]. Acta Geotechnica, 2021: 1-15.

[42] Wang T, Wautier A, Zhu J, et al. Stabilizing role of coarse grains in cohesionless overfilled binary mixtures: A DEM investigation[J]. Computers and Geotechnics, 2023, 162: 105625.

[43] Wautier, A., Bonelli, S., Nicot, F. Rattlers' contribution to granular plasticity and mechanical stability[J]. International Journal of Plasticity, 2019, 112: 172-193.

[44] Xiao P, Liu H, Stuedlein A W, et al. Effect of relative density and biocementation on cyclic response of calcareous sand[J]. Canadian Geotechnical Journal, 2019, 56(12): 1849-1862.

[45] Xu, J., Liu, J., Shi, X., Ji, F. and Zeng, Y. DEM investigation of the effect of coarse content and size ratio on the local shear behavior at gap-graded soil-structure interface. Powder Technology, 2024, 433, 119281.

[46] Xu M, Guo J. DEM study on the development of the earth pressure of granular materials subjected to lateral cyclic loading[J]. Computers and Geotechnics, 2021, 130: 103915.

[47] Yang G, Yang Q, Liu W H, et al. Critical cyclic stress ratio of reconstituted silty clay under the cyclic loading[J]. Journal of Convergence Information Technology, 2012, 7(23): 1-8.

[48] Zhao J, and Guo N. Unique critical state characteristics in granular media considering fabric anisotropy. Géotechnique, 2013, 63(8):
409 **Appendix A1: Calculation of a<sub>c</sub> and a<sub>n</sub>**

410 The fabric tensor proposed by Oda (1982) quantifies the directional distribution of contact normals.

$$
\phi_{ij} = \int_{\Omega} E(\mathbf{n}) n_i n_j d\Omega = \frac{1}{N_c} \sum_{c \in N_c} n_i n_j \quad (A1)
$$

411 where **n** is the unit vector along the normal direction of the contact plane, nᵢ and nⱼ are the ith and jth component of **n**,
412 **Ω** characterizes of the space of all the directions **n** relative to the global coordinate system, Nᵣ is the total contact
413 numbers, and E(**n**) is the distribution probability function for having a contact along direction **n**. E(**n**) can be
414 expressed as the following second-order approximation:

$$
E(\mathbf{n}) = \frac{1}{4\pi} (1 + a_{ij}^c n_i n_j) \quad (A2)
$$

415 where the second-order deviatoric anisotropy tensor $a_{ij}^c$ is symmetric and characterizes the fabric anisotropy. $a_{ij}^c$ can
416 be determined using the following expression:

$$
a_{ij}^{c} = \frac{15}{2} \phi'_{ij} \qquad (A3)
$$

417 where $\phi'_{ij}$ is the deviatoric part of $\phi_{ij}$.

418 As the mechanical anisotropy is related to normal force anisotropy, the following second-order tensor $a_{ij}^n$ is
419 adopted:

$$
a_{ij}^{n} = \frac{15}{2} \frac{\chi_{ij}'^{n}}{\bar{f}^{0}} \qquad (A4)
$$

420 where $\chi_{ij}'^n$ is the deviatoric part of $\chi_{ij}^n$, $\chi_{ij}^n$ being determined from Equation (A5). $\bar{f}^0 = \chi_{ii}^n$.

$$
\chi_{ij}^{n} = \frac{1}{N_c} \sum_{c \in N_c} \frac{f^n n_i n_j}{1 + a_{kt}^{c} n_k n_l} \quad (A5)
$$

421 Because **a**<sup>c</sup> and **a**<sup>n</sup> are deviatoric tensors by definition, it is convenient to use the following invariants to
422 quantify the degree of anisotropy:

$$
a_c = \operatorname{sign}(a_{ij}^c \sigma_{ij}') \sqrt{\frac{3 a_{ij}^c a_{ij}^c}{2}} \quad (A6)
$$

$$
a_n = \operatorname{sign}(a_{ij}^n \sigma_{ij}') \sqrt{\frac{3 a_{ij}^n a_{ij}^n}{2}} \quad (A7)
$$

423 where $\sigma'_{ij}$ is the deviatoric part of the stress tensor $\sigma_{ij}$ and sign() is the sign function. The sign function gives the
424 relative orientation of the principal direction of $a_{ij}^{\epsilon}$ and $a_{ij}^{n}$ with respect to that of the stress tensor. A positive sign

425 indicates that the major principal direction of $a_{ij}^{\epsilon}$ or $a_{ij}^{n}$ is closer to the major principal direction of the stress tensor.