To investigate the nonlinearity and energy-dissipation characteristics of particle contacts under impact, the numerical scheme employs nonlinear Hertz contact springs in both the normal and tangential branches, each arranged in parallel with a viscous dashpot, while tangential slip is governed by Coulomb friction (M’boungui et al., 2014; Horabik et al., 2017; Wu et al., 2023; Rozbroj et al., 2025). Tensile forces and contact moments are neglected; contacts carry compression only and do not impede relative rotation, consistent with the contact-kinematics schematic in Figure 2. The contact force is composed of a nonlinear elastic component and a viscous damping component, as defined in Equations 1-3, denoted as: F c = F h + F d , M c = 0 (1) F h = F n h n + F s h t (2) F d = — β n v n n — β s v s t (3)

Where F h denotes the Hertz elastic component and F d the viscous component; n and t are the unit vectors in the normal and tangential directions, respectively (Fang et al., 2025; Shi et al., 2024). v n and v s are the relative normal and tangential velocities at the contact point, and β n and β s are the normal and tangential damping coefficients. The normal branch is tensionless: F n h ≥ 0 , otherwise, the contact fails. The tangential branch follows the Coulomb criterion: | F s h | ≤ μ F n h

When the limit is reached, sliding occurs and the tangential displacement increment is decomposed as given in Equation 4: Δ δ s = Δ δ s e + Δ δ s μ (4) where Δ δ s e is the elastic part and Δ δ s μ is the sliding part.

At the contact scale, energy is partitioned into elastic strain energy E k , sliding dissipation work E μ ​, and viscous damping dissipation work E β (Mindlin and Deresiewicz, 1953; Cundall, 1988).

Elastic strain energy E _k means energy stored in the nonlinear springs,is calculated by Equation 5: E k = α h α h + 1 F n h 2 k n + 1 2 F s h 2 k s (5) where α h is a restitution-related coefficient ensuring near-constant restitution at low impact velocities (Hunt and Crossley, 1975). K n ​ and k s ​ denote the initial normal and tangential stiffnesses, which depend on the current normal and tangential force states; F n h and F s h are the Hertz normal and tangential forces, respectively.

Sliding dissipation work E μ was dissipated by tangential slip is computed by Equation 6: E μ = 1 2 F s h o + F s h n · Δ δ s μ (6) where F s h o and F s h n are the tangential forces at the beginning and end of the time step.

Viscous damping dissipation work E _β was dissipated by the dashpots is defined in Equation 7: E β = F d · δ ˙ Δ t (7) where δ ˙ is the relative translational velocity defined in the contact-kinematics analysis, and Δ t is the time-step size. The damping ratios were set to β n ​ = 0.7 and β s and kept identical for all scenarios. In this study, E β is treated as an effective rate-dependent dissipation term at the contact scale, while it may also partly serve a numerical stabilization role. Accordingly, E β is used for energy partitioning and comparative trend analyses, rather than being directly interpreted as a transferable, material-intrinsic damping property without qualification.

For the horizontal direction (Figure 6a), at 6,400 kN·m the porosities after 10 punches are 33.7% at m1, 34.8% at m2 (horizontal distance ≈0.4 m from the impact point), and 35.3% at m3 (≈0.6 m). After five squeezes with 50% gravel, m2 and m3 register 34.1% and 34.5%, respectively; with 60% gravel, m2 and m3 are 34.2% and 34.7%. At 8,000 kN·m, the corresponding porosities are 34.7%, 34.1%, and 34.2% at m2, and 35.3%, 34.5%, and 34.7% at m3 (Figure 6b). For a given energy level and gravel content, porosity along the horizontal line increases markedly with horizontal distance from the impact point. Moreover, at 8,000 kN·m the m1 porosity is consistently lower than at 6,400 kN·m, whereas m2 and m3 are broadly comparable between the two energy levels, indicating rapid attenuation of energy with horizontal distance and, consequently, significantly lower densification at mid-field and far-field locations (m2, m3) than near-field (m1).

For the 45° oblique points (m4, m5), at 6,400 kN·m the porosities after 10 punches are 32.9% and 34.3%; after five squeezes with 50% gravel they are 33.1% and 34.1%; with 60% gravel they are 33.1% and 34.2% (Figure 6c). At 8,000 kN·m, the corresponding values for m4 are 32.8%, 32.5%, and 32.6%, and for m5 are 33.9%, 33.6%, and 34.0% (Figure 6d). Compared with the horizontal direction, porosity at the 45° oblique points is lower at both energy levels, implying slower attenuation of energy along the oblique path, which is related to higher efficiency of particle collisions and energy transfer along oblique stress-wave trajectories.

For the vertical points (m6, m7), at 6,400 kN·m the porosity after 10 punches is 32.2% at m6 and 32.9% at m7. After five squeezes with 50% gravel, m6 and m7 are 32.0% and 32.1%; with 60% gravel they are 32.2% and 32.4% (Figure 6e). At 8,000 kN·m, the corresponding values are 31.5%, 31.9%, and 32.2% at m6, and 30.3%, 30.7%, and 30.9% at m7 (Figure 6f). The data show a stable decrease of porosity with depth in the vertical direction and, at the same depth, significantly lower porosity than along the horizontal and 45° oblique directions (e.g., under 8,000 kN·m with 60% gravel, m7 is far below m5 at 34.0%). This indicates the slowest attenuation of energy vertically, attributable to alignment of the vertical direction with the path of impact-force transmission, which minimizes wave-energy loss and promotes more extensive particle migration and pore compression, thereby achieving the highest degree of densification.

Under the Hertz hysteretic–viscous contact formulation, input energy during dynamic compaction is primarily transformed and dissipated through three channels—elastic strain energy, sliding work, and viscous damping dissipation (Mindlin and Deresiewicz, 1953). Using the composite energy response plots for two model input energy levels (40J and 50J) and two gravel contents (50% and 60%) (Figure 7), the post–five blow energy budgets are compared to elucidate how energy level and mixture ratio govern transmission efficiency and dissipation mechanisms.

At 40 J of model input energy, stress-wave amplitudes and particle motions remain relatively moderate. The magnitudes of E k and E β ​ are comparable, and the differences between the 50% and 60% gravel contents are small (Figures 7a,b). At the end of five blows, the slip work E μ ​ is approximately 9.5 J for the 50% mixture and 10 J for the 60% mixture, with a difference of about 0.5 J. This reflects a moderate energy input with balanced particle-scale dissipation, in which elastic compression and velocity-dependent damping jointly dominate the energy budget and variations in gravel content do not markedly alter the partitioning. When the model input energy increases to 50 J, high-energy stress waves substantially intensify particle motion and deformation. In this regime, E μ becomes more sensitive to gravel content. Within five blows, E μ ​ is about 6 J for the 50% case and about 3 J for the 60% case, indicating that poorer skeletal continuity promotes more random slip and greater dissipation through non-directional friction (Figures 7c,d). In addition, the interlocking and rolling-resistance effects associated with angular gravel are not explicitly represented by the spherical-particle idealization, which may underestimate skeleton continuity and bias the predicted sliding work toward higher values under strong impacts, particularly for the more discontinuous 50% gravel case. As shown in Figure 7c, the discontinuous skeleton at 50% gravel, coupled with reduced kinematic constraint on elastic deformation, amplifies particle vibration velocities, rendering E k in the 50% case the highest among all scenarios.

To facilitate cross-condition comparison and to clarify the engineering significance of the observed differences, the sliding work is further reported in normalized, increment-based forms. Specifically, the increment-based sliding fraction within the contact-energy budget is defined as R μ Δ ​, and the increment-based sliding ratio normalized by the external input energy on the model scale is defined as R μ i n , Δ ​, with the evaluation window taken from the end of the first blow to the end of the fifth blow. The resulting indices are summarized in Table 5. The table shows that, at the medium energy level, both mixtures exhibit comparable R μ Δ and R μ i n , Δ ​, indicating a relatively stable partitioning. At the high energy level, the 50% mixture presents a noticeably larger input-normalized sliding increment than the 60% mixture, confirming that the difference reflects a change in energy partitioning mechanisms rather than only the absolute magnitude of E μ ​. This quantitative normalization supports the interpretation that a less continuous skeleton tends to promote more frictional slip dissipation under strong impacts.

The two-dimensional μCT slices and the corresponding three-dimensional reconstruction in Figure 8 indicate that a 60% gravel content develops a more continuous load-bearing skeleton with a more regular external morphology, whereas 50% gravel exhibits localized discreteness and discontinuities that provide insufficient rigid constraint. Consequently, at high energy levels, the 50% mixture is prone to disordered, repeated slip and collision events, producing stepwise increments of E μ that are markedly larger than in the 60% case. By contrast, the 60% case shows more directional and controlled particle kinematics with less random sliding, so a larger share of the contact-energy budget contributes to compressive deformation and densification rather than friction-dominated dissipation. Energy is transmitted most efficiently along the vertical path to drive densification, while horizontal transmission decays rapidly with distance. This is consistent with the earlier “rapid-then-gradual, approaching steady” evolution of settlement with blow count and with the spatial pattern of porosity—lowest vertically and highest horizontally. In other words, vertical transmission more readily converts input into effective E μ + E k ​, whereas horizontal transmission, being more susceptible to boundary effects and heterogeneity, is preferentially transformed into E β and dissipated quickly.

Using the peak velocities from five squeezing blows marked by blue/pink boxes in Figure 10 and the corresponding directional distances, exponential engineering fits were performed for the vertical Z, the 45° oblique S, and the horizontal X directions (Figure 11), as expressed in Equations 8-10: V z z = V z 0 · e − α z · z (8) V d s = V d 0 · e − α d · s (9) V x X = V x 0 · e − α x · x (10)

The goodness-of-fit is reported using the coefficient of determination R 2 , which is annotated on each fitting curve in Figure 11. For the three principal directions, the exponential regressions provide consistently high R 2 values, supporting the suitability of the exponential attenuation form within the investigated distance range. Results indicate exponential decay in all three directions, with attenuation strength satisfying: horizontal> 45∘ oblique > vertical. This ordering is consistent with large-scale field observations at Daxing Airport, where lateral attenuation exceeds downward attenuation, indicating predominantly downward stress-wave transmission (Jifang et al., 2025). This consistency is qualitative and refers only to the directional ranking, not to matched attenuation coefficients or scaled distances.

Gravel content exerts a clear regulatory effect on the attenuation coefficients: the 60% mixture exhibits smaller vertical, 45° oblique, and horizontal attenuation, whereas the 50% mixture shows higher near-surface peaks but stronger decay in all directions (Figure 11). This behavior accords with the energy partition mechanism under the Hertz hysteretic–viscous formulation: insufficient skeletal continuity at 50% enhances random interparticle slip, elevating sliding work E μ ​ and promoting shallow and lateral dissipation, hence larger attenuation; By contrast, the 60% mixture forms more continuous force chains and a stiffer skeleton, increasing wave impedance and strengthening downward energy transmission, which manifests as slower vertical decay and suppressed horizontal diffusion. The directional attenuation is corroborated by micromechanical densification indicators: vertical porosity decreases steadily with depth and is markedly lower than that at equal distances along the horizontal and 45° oblique directions; At the same energy level, near-surface point m1 shows substantially higher porosity (lower densification) differences relative to middle and far field points m2 and m3, reflecting the rapid lateral attenuation of energy.

According to vibration stress-wave propagation theory (Singla and Gupta, 2021; Bao et al., 2025), tamping in soils can generate three representative wavefronts: (1) a quasi-one-dimensional, plane-wave-like downward propagation with limited geometric spreading within the investigated domain, where the wavefront can be approximated as locally planar; (Figure 12a); (2) a cylindrical front diffusing parallel to the ground surface (Figure 12b); and (3) a hemispherical front emanating from the pit bottom (Figure 12c). Distinct geometries induce distinct geometric-scattering terms. To simultaneously represent geometric spreading and viscoelastic dissipation of the medium, a direction-dependent exponential family is introduced: V l , θ = A 0 θ · l − n θ · ⁡ exp − α θ l + β θ l 2 (11) where A 0 θ is the directional initial amplitude, n θ ​ is the geometric-scattering exponent, and l is the path length along the wavefront in direction θ . The coefficient α θ captures first-order losses associated with viscosity and small-scale scattering, whereas β θ ​ accounts for second-order losses accumulated along the path due to multiple scattering and coupling with heterogeneity. When β θ ≈ 0 , the model degenerates to single exponential decay; When β θ >   0 , it yields an exponential with a quadratic polynomial form, enabling description of the near to far field transition from viscosity dominated to multiple scattering dominated attenuation.

From the perspective of classical stress-wave propagation in tamped soils, the three representative wavefronts reported in the literature (locally planar, cylindrical and hemispherical) mainly differ in their geometric spreading and the associated energy dilution with distance. Equation 11 is consistent with this framework by separating a geometry-controlled term l − n θ from a path-loss term exp − α θ l + β θ l 2 : the former captures the effective strength of geometric spreading for a given wavefront geometry and direction, whereas the latter summarizes irreversible losses accumulated along the path due to contact damping and interface-related scattering in the clay–gravel mixture. In this sense, the present wavefront-based formulation does not replace existing propagation theory, but provides a compact way to map the observed directional PPV decays onto representative wavefront geometries, thereby enabling direct comparison of directional ranking and gravel-content effects within a unified expression.

Importantly, this unified formulation is built upon the directional fits and ordering identified in §4.2.1: the empirically obtained α z 、 α d 、 α x ​ are regarded as directional instances of α θ . This formulation embeds the fitted directional exponential attenuation into representative three-dimensional wavefront geometries and offers a wavefront-consistent, physics-informed interpretation of the empirical trends; the parameters n θ ​, α θ and β θ ​ are treated as regression-derived effective descriptors rather than independently calibrated wavefront constants, and the formulation is therefore not presented as a fully calibrated mechanistic model.