Investigating the influence of excitation force on compaction behavior of gravel in granular blends through DEM simulation

During the construction of subgrades, the compaction of granular blends is typically achieved through vibrational techniques at varying the frequencies and excitation forces. This study investigates the impact of excitation forces, within the range of 300-600 kPa, on subgrade compaction by employing the discrete element method (DEM). Calibration tests were conducted to establish the contact parameters for the DEM models, and their reliability was verified via vibration compaction tests. The research further explores the evolution of settlement, particle motion, and contact interactions at both macroscopic and mesoscopic levels as influenced by the excitation forces. Furthermore, the influence of excitation force on the stability of skeleton frameworks was analyzed. The results indicated that the number of contacts in the final state increased consistently with the excitation force, leading to a more uniform distribution. This change contributes to a time-dependent stability within the skeleton framework, which effectively limits the movement of fine particles in the final stages and diminishes the sliding between the coarse particles. Conversely, in the initial phases, a rise in excitation force increases the stress concentration between the contacts, which increases the sliding damage to the skeleton frameworks and leads to greater compaction deformation.


Introduction
The structural integrity of cohesionless granular materials such as mixtures of sand and gravel depend on the friction and interlocking between particles [1].The mechanical properties of these materials, including strength and stiffness, are intricately linked to their compaction quality.Liu and Chen [4] demonstrated that an increase in compactness enhances the strength and stiffness of granular materials.Therefore, compactness-measured through parameters such as settlement, dry density, or porosity-is crucial for evaluating the quality of compaction [2].
The existing research has primarily conducted triaxial compression or direct shear tests to determine the influence of particle shapes on the shear strength [4] and resilient deformation [5] of granular blends.These methodologies predominantly assess the overall effect of particle shape on mechanical properties at a macroscopic level, which were constrained by experimental conditions.To conduct a thorough analysis that includes both macroscopic and mesoscopic mechanical characteristics [6], the discrete element method (DEM) has been recognized as an effective tool.This method is extensively used for investigating granular materials, providing detailed insights into the influence of particle shapes on mesoscale properties [7].1332 (2024) 012021 IOP Publishing doi:10.1088/1755-1315/1332/1/012021 2 Numerical investigations into particle shape effects often involve generating polyhedral models representing various typical particle shapes [8] or utilizing advanced three-dimensional (3D) reconstruction techniques to mimic the actual particle geometries [9].Lai et al. [10] assessed how 3D particle shapes influence the compaction properties of granular blends under vibrational loads.Initial examinations of the impact of vibration parameters on the packing density of regular and irregular particles under vibrational loads have been conducted by [11], yet a detailed quantification of these effects remains unaddressed.These investigations primarily used simplified particle shapes, neglecting the natural irregularities of real particles.Thus, a quantitative study on how excitation force affects the compaction behavior of granular blends from a mesoscale perspective is essential to bridge these gaps in knowledge.

Statistical analysis of gravel shapes
For the accurate modeling of gravel blends that reflect the authentic shapes of individual gravel particles, an analysis of the shape characteristics of these particles was undertaken.The particle shape is evaluated at various geometric levels.Prior research has identified three key criteria for quantifying particle shape: overall shape, roundness, and surface texture [12,13].Considering that the surface texture is effectively represented by the friction coefficient value in DEM simulations, and this study aims to capture the wide variety of particle shapes, the overall shape emerges as the most suitable criterion for characterization.Regarding the general shape configuration, shape indices are calculated by comparing the dimensions of the three principal axes.Elongation (EI) is measured as the ratio of the secondary axis length to the major axis length, whereas Flatness (FI) is the ratio of the minor axis length to the secondary axis length.The EI and FI values range from 0 to 1.A value nearing 0 indicates a more elongated or flattened particle, whereas a value approaching 1 suggests a particle with nearly identical diameters.
In this research, 230 gravel particles were initially converted into point cloud models using a 3D laser scanner, and their shape characteristics were quantified.The statistical probability distributions of EI and FI were examined, as depicted in Figure 1, displaying a broader range of EI values compared to FI values.Consequently, EI was selected as the primary descriptor for selecting gravel templates to produce both coarse and fine particles within the model.Based on the probability statistics of EI, gravel particles were created.Specifically, six gravel particles were selected as templates, each with a consistent FI value of 0.75 but with varying EI values of 0.48, 0.58, 0.64, 0.73, 0.81, and 0.92.

Contact parameter calibration
In this discrete element simulation, the Hertz-Mindlin (No Slip) model was employed, focusing on primary contact parameters such as Young's modulus, coefficients of restitution, and coefficients of friction.The accurate calibration of these parameters is essential for the validity of the results [14].To calibrate the microscopic parameters, compression, collision, and inclined plate tests were conducted by performing each test thrice per group and utilizing the average values.The calibration of Young's modulus (E) required adjusting its value in the DEM model to align with the outcomes of physical tests.This calibration process involved a spherical particle in both physical and DEM tests, ultimately setting E to 4.9 × 10  after several iterations.An inclined plate test was conducted [15] to calibrate the particle-wall friction coefficient ( ), resulting in an average angle ( ) = 24.2°and ( ) of 0.45.The particle-particle friction coefficient ( ) was calibrated similarly, with an average angle ( ) = 28.4°,resulting in ( ) of 0.54.The restitution coefficients ( ) and ( ) were adjusted using single-particle and double-particle collision tests [15].For ( ), a value of 0.65 was determined from a single-particle collision test, whereas for ( ), a value of 0.47 was determined from a double-particle collision test, as depicted in Figure 2. The calibrated contact parameters and their values used in the subsequent DEM models are summarized in Table 1 based on the calibration process and material manuals.

Simulation conditions for DEM models
This study utilized six gravel particles with diverse elongation (EI) values to generate both coarse and fine particles.For the creation of fine particles, the particle with an EI value of 0.92 was selected because of its computational efficiency.The generation of effective particles was influenced by factors such as the smoothing value ( ), minimum sphere radius ( ), and software algorithm parameters.A equilibrium between computational efficiency and accuracy was achieved with ( ) = 5 and ( ) = 0 mm. Figure 3 illustrates the number of sub-spheres in the six gravel templates.To optimize the computational cost and timestep of the DEM model while maintaining reasonable limits, the minimum particle size for numerical analysis was set at 5 mm, corresponding to 1/12 th of the maximum particle size [16].Consequently, particles with a dimension of 4 mm were categorized as fine particles.The simulation initiated with the random generation of fine and coarse multi-sphere particles within a steel-walled cylindrical container ( ×  = 155 mm × 190 mm) , as displayed in Figure 3. Subsequently, the particles underwent a pouring process, settling under the influence of gravity.Interactions and collisions with the container walls dissipated energy, reducing the velocity of particle movement.A stable packing configuration was achieved once all particles attained equilibrium positions, typically within approximately 0.5 s.At 1.0 s, vibration was introduced to form an initial packing structure [17].
To apply pressure to the granular blend, a load plate ( × ℎ = 152 mm × 10 mm) is used in two stages: preloading at a constant 100 kPa and vibration loading with a sinusoidal waveform  = − − + sin   + in the Z-direction within the container.The vibration period lasted for 10 s, succeeded by 1.0 s of stabilization.The boundary effects from the specimen and particle sizes are considered negligible when the ratio between the specimen size and particle size was greater than 5 [18].In this study, the steel-walled cylinder was approximately four times larger than the largest particles, which minimized the boundary effects.

DEM model validation
Utilizing an average force of 450 kPa as a reference, the validation process applied forces of 400 kPa and 500 kPa.Through detailed comparisons between numerical simulations and corresponding physical tests utilizing these forces, significant insights were obtained.Both simulations and measurements from the multifunctional testing machine, as illustrated in Figure 4, revealed a consistent pattern.Initially, there was a rapid increase in settlements, which then tended to stabilize.Notably, settlements observed during physical testing rose and stabilized more quickly than those in the simulations.This difference was attributed to the greater diversity in particle shapes in the physical tests, leading to enhanced interlocking among particles and a faster formation of a stable framework.The evolution curves of both the tested and simulated settlements are plotted in Figure 5 (b) along with the relative error calculated as the ratio of the difference between tested and simulated settlements to the tested settlements.Despite certain differences between the test and simulation outcomes, these were primarily prominent during the strong deformation phase, with negligible effects observed during the stable deformation phase.Therefore, the DEM model accurately reflects the actual deformation response under vibrational loads, albeit with minor discrepancies during specific deformation stages.

Evolution Permanent deformation
The variations in settlement caused by different excitation forces serve as indicators of the progression of permanent deformation.Initially, all samples subjected to testing exhibited a rapid increase in settlement, which then gradually increased before stabilizing after two seconds, as depicted in Figure 5 (a).To further investigate the effect of force on settlement evolution, Figure 5 (c) presents the initial settlement during compaction, the final settlement post-compaction, and the settlement reduction.The data indicate that both initial and final settlements, as well as settlement reductions, decrease and then increase with higher forces, suggesting that lower forces result in less compacted samples.The difference between the initial and final settlements-termed as settlement reduction-is detailed in Figure 5 (c), highlighting that settlement reduction escalates with increasing forces, with the minimal point recorded at 400 kPa.

Evolution of particle motion
Compaction deformation occurs when particles move through translation and rotation.Figure 6 (a) and Figure 6 (b) depict the translational and rotational velocities of coarse and fine particles throughout the loading process.As the loading time increased, the translational velocity decreased rapidly and then tends to stabilize after approximately two seconds, indicating an inverse correlation with settlement changes.Conversely, the decline in rotational velocity of coarse particles was comparatively modest [10].Both translational and rotational velocities demonstrated a strong inverse relationship with settlement, suggesting that compaction deformation involves both translation and rotation of particles.
Figure 6 (c) and Figure 6 (d) display the variations in average translational and rotational velocities as influenced by force.The rotational movement of fine particles consistently exceeded that of coarse particles, whereas the translational movement of coarse particles generally exceeded that of fine particles.Fine particles play a critical role in filling the voids between coarse particles, which aid in the formation of densely packed samples.With an increasing force, both coarse and fine particles exhibit an uptick in velocity, although the rotational movement of coarse particles remains largely unaffected.

Mesoscopic contact evolution and analysis
This section examines into the impact of force on contact characteristics, including the average number of contacts and the spatial distribution of contact forces.

Contact numbers
The number of contacts is a crucial metric for quantifying the internal structural characteristics of granular materials, traditionally defined as the average number of particles in contact with a given particle.The average contact number (Z) denotes the mean number of contacts per particle [19].The average contact number can be numerically calculated as follows: where  ,  , and  denote the numbers of  ,  , and  , respectively. and  represent the numbers of fine and coarse particles, respectively.Figure 7 (a) illustrates the evolution of Z over time during the loading phase.Initially, Z rapidly increases, then slows and stabilizes after two seconds, reflecting the pattern of settlement evolution.Figure 7 (b) presents variations in Z under different forces, outlining the initial, final, and reduction stages.As force increases, the average contact numbers in the final state consistently rise, whereas those in the initial state initially decrease then increase.The reduction stage shows growth and is consistent with the trend observed in settlement reduction.
The contacts were categorized into three types based on the particle interaction: coarse-coarse contacts (CC contacts), coarse-fine contacts (CF contacts), and fine-fine contacts (FF contacts), denoted as ZCC, ZCF, and ZFF, respectively.Figure 7 (c) and Figure 7 (d) illustrates the variations in the number of these partial contacts over time under loading conditions.All three types of partial contacts experienced a rapid increase followed by gradual stabilization after 2 s, reflecting a consistent pattern with settlement development.This pattern indicates that contact relationships tend to stabilize as the sample compacts.During the loading phase, ZCC, ZCF, and ZFF numbers initially increased and then decreased with increasing force, suggesting a significant influence of excitation force on CC, CF, and FF contacts.The hierarchy of contact numbers among these contacts is observed as ZCF > ZFF > ZCC.
Figure 7 (e) and Figure 7 (f) display the fluctuations in partial contact numbers with varying force.The initial states of ZCF and ZFF decreased, whereas ZCC consistently increase with the intensifying force.Furthermore, the final states of ZCC and ZCF numbers consistently increased, whereas the ZFF numbers decreased monotonously.The contact number reduction depicted in Figure 7 (f) shows that increases in ZCC, ZCF, and ZFF correlate with increases in settlement reduction, indicating a correlation between all contact types and settlement variation.

Contact force spatial distribution
The magnitude of contact forces serves as a mesoscopic metric to assess load resistance, which reflects the difficulty of compaction.In this section, we analyze the contact forces of the three partial contacts and the distribution of the sliding and non-sliding contacts.The variations in contact forces for the three types of partial contacts are presented in Figure 8, including the average contact forces for CC contacts (FCC), CF contacts (FCF), and FF contacts (FFF).The order of contact forces, FCC > FCF > FFF, suggests that the primary strength originates from CC contacts.Additionally, FCF and FFF demonstrate a descending order in response to excitation force, indicating that excitation force primarily affects the contact forces involving coarse particles with minimal impact on forces between fine particles.
The sliding contacts are governed by Coulomb's friction law.The sliding ratio  can be defined as  = | | , where  and  denote the shear and normal forces, respectively. denotes the friction coefficient between particles.A sliding contact is considered to occur when  ≥ 0.9999 [9].The percentage of sliding (SCP) and be expressed as follows: where  denotes the number of sliding contacts.
In Figure 9 (b), a negative linear correlation is observed between the SCP and settlement.Notably, the SCP decreases as settlement increases, signifying that the granular material becomes more resistant to shear forces over time.This increased resistance could be attributed to the enhanced anisotropy of contact forces, which impedes the sliding of particles past one another.
The distribution of the sliding contacts within the X-Z plane is illustrated in Figure 10, along with the SCP for each sample.Interestingly, the SCP decreased as the loading time extends.During the initial loading times of 2.5 s and 5.0 s, the SCP increased with the excitation forces, attributed to the heightened collision frequency and intensity among particles, leading to an increase in relative sliding.Conversely, at loading times of 7.5 s and 10 s, SCP decreases under stronger excitation forces as the particles are forced into denser packing arrangements, limiting the space available for sliding.

Discussion
Further exploration into the impact of excitation forces on the fluidity and structural stability of the skeletal framework involves analyzing the interplay between particle movement and contact.Fine particles demonstrate significantly greater rotational motion compared to coarse particles, whereas coarse particles exhibit superior translational motion.Moreover, ZCF and ZFF are more profoundly impacted than ZCC, highlighting that changes in particle movement predominantly result in alterations in contact points, especially among contacts involving fine particles.Additionally, the sequence of contact force interactions is stated as follows: FCC > FCF > FFF, suggesting that coarse particles primarily bear the brunt of transmitting contact forces.Consequently, the skeletal framework is primarily composed of coarse particles, with fine particles integrated within it.The variation of the SCP with excitation force is depicted in Figure 9 (a).SCP exhibits a fluctuationinitially increasing and then decreasing-as excitation force intensifies, suggesting a time-dependent behavior of granular particles.Initially, particles are more susceptible to sliding as the forces increase, but over time, particularly at 7.5 and 10 s, they begin to compact, possibly because of particle rearrangement, settling, or other gradual compaction processes.The initial instability can be ascribed to the uneven distribution of contacts that emerges with the increasing excitation forces, where gravel particles located at the stress accumulation points tend to slip more readily, disturbing the stability of the skeletal frameworks.Thus, alterations in excitation forces significantly influence the stability of these frameworks by heightening stress concentrations at contact points.

Conclusions
The objective of this study was to explore the impact of excitation force on the compaction characteristics of gravel in granular blends, employing real-shaped gravel with a specific particle size distribution and creating DEM models.The study produced several key findings:  The initial, final, and reduction settlement increase with the compaction force, which can be attributed to the periodic motion of the gravel particles that influences the interaction between the coarse and fine particles and affects the compaction capability of the granular blends. The skeletal framework is primarily composed of coarse particles, with fine particles occupying the interstitial spaces.An increase in force results in the framework exhibiting time-dependent stability.In its stable phase, the framework restricts the movement of fine particles, leading to a more uniform and stable granular arrangement. Variations in excitation forces affect the stability of the skeletal frameworks by enhancing stress concentrations at the contact points.This effect is observed alongside the time-dependent stability of the material during compaction, where stress concentrations become more pronounced either at the beginning or end of the compaction process.Simplifying the PSD curve could mitigate model compressibility.Future research should strive to develop an efficient DEM model that encompasses complete PSD curves.Although flexible boundaries may facilitate more realistic transverse deformation, this study does not conclusively determine if the effects of rigid boundaries were entirely neutralized, due to the limited boundary size.

Figure 2 .
Figure 2. Schematics of discrete element simulation calibration test

Figure 3 .
Figure 3. Formation of the blend and the application of vibrating loads

Figure 4 .
Figure 4. Testing machine and tested material

Figure 5 .
Figure 5. Settlement evolution with time and force: (a) Evolution of settlement deformation, (b) evolution curves of test and simulation, (c) evolution curves of initial, final and reduction settlement.

Figure 6 .
Figure 6.Velocity evolution with time and force: (a) Variations in translational velocity over time as loading progresses, (b) variations in rotational velocity over time as loading progresses, (c) variations in translational velocities with force, (d) variations in rotational velocities with force.

Figure 7 .
Figure 7. Contact number evolution with time and force: (a) Variations in overall contact number as loading time progresses, (b) variations in initial, final and reduction contact number with force, (c) and (d) Variations in partial contact numbers as loading time progresses, (e) and (f) Variations in partial contact numbers and reduction contact with force.

Figure 8 .Figure 9 .Figure 10 .
Figure 8. Evolution curves of partial contact forces: (a) CC contact forces (b) CF and FF contacts

Table 1 .
Contact parameters and material parameters in DEM models.