Numerical study on phase difference-driven mechanical behaviors of gravel–structure interface under three-dimensional circular cyclic shear

Abstract

Soil–structural interfaces, such as those between foundations, tunnels, retaining walls, and slopes, serve as critical zones for stress and deformation transfer in geotechnical systems. Although extensive research exists on two-dimensional interface mechanics, the zonal characteristics, non-coaxial behavior, and reversible volumetric deformation mechanisms under three-dimensional (3D) loading remain poorly understood. To address this gap, a series of discrete element method (DEM) simulations focused on the phase-dependent response of shear stress, normal displacement, and microstructural interactions were conducted in this study. The influence of the tangential displacement phase difference D_p on the 3D cyclic shear behavior of gravel–structure interfaces under circular shear paths was investigated based on DEM. The controlling role of shear path geometry in governing non-coaxial behavior and reversible volumetric deformation is explicitly revealed. Specifically, it is demonstrated that stable, hysteretic stress–displacement loops with minimal volumetric hysteresis are promoted by circular shear paths (D_p = π/2, 3π/2), whereas abrupt stress reversals and amplified dilatancy-recompression cycles are induced by bidirectional linear paths (D_p = π, 2π). Furthermore, periodic oscillations of non-coaxial angles (α₃) are observed, transitioning from 90° under pure rotation to peak fluctuations exceeding 30° under phase-shifted shear. These mechanistic understandings are highlighted as being crucial for the advancement of 3D interface constitutive models.

Graphical abstract

Keywords

gravel–structure interface / cyclic shear / non-coaxial behavior / reversible volumetric deformation / phase difference / microscopic responses

1 Introduction

In the field of geotechnical engineering, solid structural materials such as foundations, tunnels, retaining walls, and slopes all exhibit contact interfaces with soil [1–4]. As the critical medium for transmitting stress and deformation between soil and structures, these interfaces play an irreplaceable role in the analysis and design of soil–structure interaction systems. Under strong confinement conditions imposed by structures, the shear characteristics of geotechnical materials often dominate the strength and deformation behavior of composite systems [5]. The direct shear test is a commonly used method in geotechnical engineering to evaluate the shear strength of soil, effectively simulating simple shear modes and plane strain conditions frequently encountered in engineering scenarios such as slopes and embankments. By obtaining the shear stress–displacement relationship, this test can reveal the mechanical characteristics of shear behavior at the macroscopic level.

Currently, numerous scholars have conducted systematic laboratory tests to investigate the monotonic and cyclic characteristics of soil–structure interfaces, covering sand–structure interfaces [6–9], clay–structure interfaces [10–12], and gravel–structure interfaces [13–15]. Most indoor interface tests adopt a two-dimensional (2D) loading mode to explore the 2D mechanical behavior of sand–structure interfaces [16–18]; however, research on 2D tests for clay–structure [10,11] and gravel–structure [19] interfaces remain relatively limited. Some researchers have further performed three-dimensional (3D) behavior tests on sand–structure and gravel–structure interfaces using customized testing devices, focusing on normal and shear responses in three orthogonal directions [13–15]. Existing studies indicate that the key controlling parameters influencing the shear resistance of structure–soil interfaces include structural surface roughness (affected by the morphology, number, and height of surface grooves [20]) and the physical-mechanical properties of soil (e.g., particle size distribution [21], water content [22], cohesion [7], internal friction angle [7], and porosity [15]).

Although direct shear tests can intuitively reveal the mechanical response characteristics of soil–structure interfaces under different working conditions, a deeper understanding of the essence of shear behavior still requires an investigation into the microscopic interaction mechanisms between particles. The discrete element method (DEM) exhibits unique advantages in exploring the meso-mechanical mechanisms of shear tests for granular soils, effectively circumventing the limitations of traditional laboratory tests, such as long sample preparation time, high cost, and restricted vertical stress range [23]. A core strength of DEM lies in its ability to accurately track the motion trajectories of particles near the interface, providing a basis for analyzing the evolution of soil structure at the meso-scale. This method supports refined meso-scale analysis, enabling systematic capture of force network topological characteristics and contact mechanical behaviors between particles––features that are difficult to observe in traditional laboratory tests [24]. Compared to conventional tests, DEM allows cost-effective simulation of different sample sizes and gradations, effectively mitigating the high costs and time-consuming nature of laboratory sample preparation. Furthermore, its flexibility in adjusting stress paths and boundary conditions provides a more comprehensive mechanistic interpretation of soil behavior under shear loading.

Numerous scholars have utilized the DEM to simulate the interfaces between granular materials and structures, providing effective means to reveal the internal action mechanisms. For instance, Bernhardt et al. [25] employed the Hertz–Mindlin contact model for 3D DEM modeling and found that the results of direct shear tests are significantly affected by the ratio of sample size to particle number. They proposed that a minimum of 10000 particles should be included in the sample to ensure simulation accuracy. Li et al. [22] constructed 3D numerical models using DEM to investigate the coupling effects of particle shape and interface roughness on shear behavior. Their findings revealed that the interaction between particle shape and interface roughness significantly influences the shear strength of particles, with a critical threshold identified. Nina and Athanasopoulos-Zekkos [26] conducted DEM simulations on a direct shear model and validated the simulation results against steel ball direct shear tests. They emphasized the importance of contact constraints between fixed steel balls and planar boundaries to prevent particle sliding and rolling, thereby enhancing the model’s authenticity. Zeraati-Shamsabadi and Sadrekarimi [24] utilized the DEM to simulate specimens with varying diameters under identical porosity conditions. By applying different scaling factors to the particle gradation, they generated diverse particle size distributions, thereby identifying the combined influences of the median particle size (D₅₀) and the diameter-to-height ratio (D/H). Based on this analysis, they proposed the concept of limiting the D₅₀ and D/H ratio to mitigate boundary effects in specimens. Zhao et al. [27] investigated the influence of particle shape on the shear strength of granular materials, the evolution of inter-particle contact networks, and force transfer mechanisms through DEM simulations.

In summary, existing studies have established the critical influences of specimen size, particle shape [28], and boundary conditions on the mechanical response of soil–structure interfaces. However, limited by the constraints of testing equipment, prior research has primarily focused on the shear characteristics of interfaces under 2D loading conditions [29–31]. In practical engineering scenarios, however, soil–structure interfaces are often subjected to 3D loading, involving the coupling of normal constraints in two orthogonal tangential directions and synchronous shear loading [32–34]. 2D interfaces, by nature, represent simplified representations of real interfaces and cannot fully capture the core characteristics of the interface’s 3D mechanical behavior (e.g., non-coaxial shear deformation behavior). Although early studies by relevant scholars have conducted a series of 3D tests to preliminarily investigate the 3D mechanical behavior of sand–structure interfaces, their research has not deeply investigated key issues such as the zonal characteristics, action mechanisms, and non-coaxial effects of normal volumetric deformation [35]. Feng and Zhang [13,15] and Rehman and Zhang [35] and Liu et al. [36] have studied the influence of shear paths on the 3D mechanical properties and normal boundary conditions of coarse-grained soil–structure contact interfaces, but the influence mechanism of complex shear path variations (e.g., path nonlinearity, multi-directional shear coupling) caused by different tangential shear displacement phase differences under circular shear paths on the mechanical behavior of interfaces remains unreported. Such evolutionary laws of complex shear paths hold significant engineering importance for the accurate prediction of interface mechanical behavior, rational determination of structural design parameters, and safety and stability assessment.

Although laboratory tests provide important references for studying the mechanical properties of soil–structure interfaces, numerical simulation techniques––particularly the DEM––offer a cost-effective and refined research tool for investigating these influencing factors [37]. Notably, the extensive application of gravel in large-scale engineering projects due to its excellent mechanical properties and wide engineering adaptability, systematic research on the complex shear paths of gravel–structure interfaces under 3D loading conditions remains relatively limited. While some studies have considered multi-directional loading, a systematic investigation into the fundamental control exerted by the tangential displacement phase difference on the microscopic and macroscopic response is still lacking. Furthermore, to reveal the underlying mechanisms of interfacial shear behavior, the influence laws of complex shear paths on the meso-mechanical responses of interfacial shearing require further exploration.

Against this backdrop, a series of direct shear tests on gravel–structure panel interfaces were conducted based on DEM in this study, with a focus on investigating the influence laws and meso-scale action mechanisms of changes in tangential shear displacement phase difference on the mechanical behavior of gravel–structure interfaces under circular shear path conditions. A DEM numerical model for 3D gravel–structure interfaces is first developed and validated. The influence characteristics of complex shear path evolution induced by different tangential shear displacement phase differences on the interfacial mechanical behavior were then analyzed, with particular attention given to discussing the tangential stress–displacement response laws, reversible normal volumetric deformation characteristics, and non-coaxial phenomena of the interface. The inherent mechanisms of interfacial mechanical responses were finally uncovered from a meso-scale perspective, especially the meso-scale characterization of contact weak surfaces under 3D loading conditions. Key support for the theoretical development of 3D shear mechanical responses of gravel–structure interfaces and engineering safety evaluation was provided by this work.

2 Numerical simulation schemes

Table 1 summarize the DEM simulations carried out in this study. A total of 19 simulations were conducted with varying normal stress conditions and phase differences. To validate the applicability of mesoscopic contact parameters and the accuracy of the proposed model, simulations of direct shear test along circular shear paths at the gravel–structure interface under varying normal stress conditions (Tests 1–4) were performed in this study. Additionally, by establishing circular shear path configurations and accounting for phase differences in tangential shear displacement variations, the influence of different shear paths (Tests 2, 5–19) on the 3D mechanical properties of the gravel–structure interface were studied, thereby providing a theoretical basis for practical engineering safety assessment.

The computational implementation of the circular shear path is illustrated in Fig. 1(a), where the phase difference between tangential displacements is set to π/2. For the circular cyclic shear tests, the interface must first undergo initial shearing along the linear path ① shown in Fig. 1(b) until reaching the target shear displacement amplitude (u_my). Subsequently, cyclic shearing is performed following the circular trajectory ②-③-④-⑤-② in repeated cycles.

In the investigation of the phase difference D_p in tangential shear displacement, a constant normal load condition was maintained at the interface, with a uniform normal stress of σ_n = 400 kPa applied across all experimental trials. The tangential displacement phase difference D_p is defined as π/2 + δ_φ, δ_φ represents incremental phase shifts (n∙π/m) between orthogonal shear components, which imply different phase offsets, where n and m are integers that define the specific magnitude of the phase offset between orthogonal shear components in the context of 3D cyclic shear tests on gravel-structure interfaces. The phase difference of tangential shear displacement was varied according to the predefined circular path parameters configured in the computational setup, and the shear path is set as shown in Eq. (1). For the shear testing protocol, the interface was first subjected to an initial shearing process along the linear trajectory depicted in Fig. 2, achieving the target shear displacement amplitude (u_my). Subsequently, cyclic shearing was performed along the curved path under the same normal stress conditions. Obviously, when D_p = π/2 + π/2 and π/2 + 3π/2, the shear path of the structural plane degenerates into two-way linear shear [15].

(1a)\begin{array}{c}{\mathit{u}}_{\mathit{x}}=2\text{π}\mathit{R}\cdot \sin (2\text{π}\mathit{t}),\end{array}

(1b)\begin{array}{c}{\mathit{u}}_{\mathit{y}}=2\text{π}\mathit{R}\cdot \sin (2\text{π}t+D_{\mathrm{p}})=2\text{π}\mathit{R}\cdot \cos (2\text{π}t+{\delta }_{\text{φ}}).\end{array}

3 Discrete element modeling for three-dimensional gravel–structure interface

DEM has emerged as a highly effective numerical simulation technique for modeling discrete materials, particularly in accurately reproducing the shear behavior of granular soils. This computational approach is extensively utilized for simulating element-scale tests under controlled laboratory conditions in soil mechanics’ research. In the present study, to allow for a clearer representation of the shear behavior and provide valuable insights into the granular material’s response under shearing conditions, numerical simulations of a direct shear test were conducted using the PFC3D 5.00 software package. The system of nonlinear equations governing particle interactions was solved through an explicit dynamic analysis scheme. Within this scheme, the equations of motion for each discrete particle were solved using a central-difference time integration method. Subsequently, incremental displacements were calculated, and contact forces were updated at each time step based on the defined force–displacement laws.

3.1 Generation of discrete element method models

The DEM models were constructed within a cylindrical mold featuring rigid top and bottom platens, replicating the boundary conditions characteristic of conventional laboratory direct shear tests for gravel–structure interfaces, as shown in Fig. 3. The structural panel is positioned atop the gravel specimen to mitigate container wall friction effects and ensure that the measured normal stress exclusively represents the interfacial stress. The cylindrical mold was designed with frictionless properties, whereas the bottom platen was assigned a significantly high friction coefficient to ensure no slippage occurred at the platen-particle interfaces. Similarly, the surface roughness of the structural panel was numerically simulated by modifying the friction coefficient parameter on the upper platen of the numerical model. This approach replicates the methodology employed in laboratory direct shear tests, where surface irregularities are typically introduced through ridges, fins, or protruding pins incorporated into both the upper and lower platens.

The homogeneous gravel specimen comprises particles with diameters ranging from D_min = 5 mm to D_max = 16 mm (where D_min and D_max denote the minimum and maximum particle diameters, respectively), exhibiting a mean particle diameter of D₅₀ = 9.0 mm (specific gravity of soil solids G_s = 2.65), the particle size distribution of the sample is shown in Fig. 4. Based on a porosity of 0.35, particles were spatially distributed in a random arrangement, as depicted in Fig. 3(a). The specimen dimensions, with a diameter of 500 mm and a height of 250 mm, exceed 10D_max of the gravel to ensure boundary effects are minimized. The structural panel has a slightly larger surface area compared to the gravel specimen, ensuring that the contact area at the interface remains consistent during the entire shearing process. Consequently, the normal displacement can be used to accurately quantify the volumetric deformation occurring at the interface. The micro-mechanical parameters of each part in the numerical simulation are shown in Table 2.

The models were loaded with the target vertical stress (σ_n) through controlled displacement of the direct shear test mold’s top platen, achieved via a wall-servo mechanism, as shown in Fig. 3(b). The effectiveness of this pre-equilibration stage was evaluated by tracking both the specimen’s vertical displacement and the interparticle unbalanced force ratio. Following compression, the consolidation void ratio (e_c) of the specimen was determined based on its measured vertical displacement. These procedures guaranteed proper sample equilibration prior to initiating the shearing phase. During the shearing process, normal pressure is applied to the structural panel through the normal loading system, while shear deformation is induced by displacing the shear box via the horizontal loading system.

3.2 Parameter calibration

Upon comprehensive evaluation of the linear model and the linear contact bond model, the rolling resistance linear model was ultimately selected for interface simulation [37]. This model demonstrates superior capability in accurately representing the contact behavior between geological materials. Specifically, it was employed to simulate particle interactions in proximity to the gravel–structure interface. The rolling resistance linear model incorporates a rolling resistance mechanism into the conventional linear model framework, effectively addressing energy dissipation and rotational resistance during particle movement, as shown in Fig. 5. This innovative approach enables the accurate reproduction of constitutive behavior for irregular soil particles through the use of simplified circular particles. Consequently, this methodology eliminates the need for repetitive evaluation of contacts between complex-shaped particles, significantly enhancing computational efficiency while maintaining simulation accuracy.

For these purposes, the rolling resistance linear model is employed to calibrate the contact model parameters through systematic comparison with experimental data [15]. The shear simulations were performed at different normal pressures (σ_n = 200, 400, 700, 1000 kPa), each specimen was sheared up to a shear displacement (u) of 40 mm. The test procedure involves sequentially applying normal loads as concentrated forces at the interface of the rigid gravel structure, followed by controlling the horizontal displacement of the shear box at a set shear rate to induce shear loads, as shown in Eq. (2). A global damping coefficient of 0.7 was implemented to accurately model energy dissipation, thereby effectively suppressing nonphysical oscillations of fine particles during DEM simulations. Throughout the shearing process, critical parameters including interfacial pressure, structural panel displacement, and shear box movement were continuously monitored to characterize both the shear stress–strain behavior and volumetric strain–shear strain response of the tested materials.

(2a)\begin{array}{c}{\gamma }_x=2\text{π}R\cdot \cos (2\text{π}t),\end{array}

(2b)\begin{array}{c}{\gamma }_y=2\text{π}R\cdot \cos (2\text{π}t+D_{\mathrm{p}})=-2\text{π}R\cdot \sin (2\text{π}t+{\delta }_{\text{φ}}).\end{array}

The direct shear test of the structural interface was numerically simulated by varying relevant parameters, ultimately obtaining an optimal parameter set (Table 3) that demonstrates excellent agreement with experimental results. The best fitting curves of shear stress–displacement response obtained by DEM compared with those obtained in the laboratory [15] were presented in Fig. 6. The numerical simulation results of the gravel–structure interface under varying normal pressures reveal that the shear stress initially increases progressively before reaching its peak value, after which it stabilizes with further tangential displacement. Notably, minimal strain-softening behavior is observed. The quantitative agreement was excellent, with a coefficient of determination (R²) of 0.999 and a root-mean-square error (RMSE) of 6.37 kPa for the peak shear stress across all normal stress levels. Furthermore, the peak stress ratio under different normal stress conditions remains consistent at approximately 0.65, aligning well with laboratory test results. This correlation confirms that the selected contact parameters effectively capture the mechanical interactions between the interface and soil particles.

3.3 Model validation

Prior to implementing the developed numerical model and calibrated contact parameters in the current DEM analysis, a series of interface shear tests (Tests 1–4) under varying normal stress conditions (σ_n = 200, 400, 700, 1000 kPa) along circular shear paths at the gravel–structure interface were performed in this study. These numerical investigations were conducted to validate the DEM model’s capability in accurately simulating the macroscopic mechanical response of the interface system. The calibrated model parameters were subsequently implemented in the circular path shear simulation process.

The shear stress–displacement responses (τ_x–u_x, τ_y–u_y) of the circular path shear of the interface under different normal stress conditions were shown in Fig. 7(a). The simulation results of the 1st, 10th, and 50th cycles are selected for display. It can be seen that except for the first cycle which is rather special, the tangential stress–strain relationship curves of the remaining cycles basically overlap. When the normal stress is 400 kPa, the peak tangential stress of the contact surface can reach about 260 kPa. Figure 7(b) illustrates the evolution pattern of normal displacement accompanying cyclic shear deformation at the interface under varying normal stress conditions. The experimental results reveal that under triaxial shear conditions, the cumulative shear displacement at the interface exhibits a monotonically increasing trend characterized by a progressively diminishing strain rate. This displacement progression indicates a transition from initial shear localization to distributed shear deformation regimes, with the shear strain rate demonstrating a gradual attenuation as cyclic loading progress. Notably, under a normal stress of 400 kPa, the interface exhibits a total shear displacement accumulation of 21.43 mm after 50 cyclic loading iterations, demonstrating a stabilized strain response typical of dense granular materials under prolonged shear loading. The results maintain a good consistency with the indoor shear test results of Feng and Zhang [15].

The experimental shear response data presented in Fig. 8 demonstrates the cyclic shear behavior of the interface under varying normal stress conditions through three representative loading cycles (1st, 10th, and 50th). The experimental results reveal that the peak tangential stress maintains a stable proportionality of approximately 0.65 relative to the applied normal stress across all investigated normal stress levels. Furthermore, the cyclic shear stress–displacement hysteresis loops exhibit consistent behavioral patterns, characterized by identical evolutionary trajectories in both shear stress magnitude and displacement accumulation irrespective of cyclic loading frequency. This cyclic stability aligns with the shear behavior observed under cruciform shear path conditions as previously documented [20]. Notably, quantitative analysis indicates a direct correlation between normal stress magnitude and interfacial shear stiffness. The initial shear stiffness (quantified as the secant modulus during the first loading branch of each cycle) demonstrates a positive correlation with increasing normal stress, exhibiting a 15%–22% stiffness enhancement per MPa normal stress increment. This stress-dependent stiffening response corroborates the interface hardening mechanism proposed by Feng et al. [38], while simultaneously demonstrating superior cyclic stability compared to rate-dependent shear models.

The 3D DEM numerical model for gravel–structure interface characterization was developed and validated through systematic simulation. Essentially, the macroscopic mechanical behavior under specified confining pressure conditions was demonstrated to exhibit high predictive accuracy. The validity of both the proposed DEM framework and the calibrated contact parameters was substantiated by the modeling results.

4 Results and discussion

The influence of phase angle variation on 3D cyclic shear behavior at gravel–structure interface was investigated through systematic parametric studies. Controlled boundary conditions were established to maintain constant shear displacement in the principal X-direction, while incremental phase shifts (δ_φ) between orthogonal shear components were progressively introduced. A total of 16 cyclic shear simulations (Tests 2, 5–19) were executed under these parameterized conditions, with δ_φ varied from π/6 to 11π/6 radians. This comprehensive parameterization allowed for examination of shear response evolution across extended phase difference regimes through controlled experimental protocols. The simulation of direct shear test was conducted on the specimen interface subjected to a maintained normal stress of 400 kPa until achieving a horizontal displacement increment of 20 mm, measured as the differential displacement between structure panel and shear box (u_x or u_y).

4.1 Shear stress–displacement hysteretic response

To investigate the impact of tangential displacement phase lag on the circular shear behavior of interfaces, a comprehensive analysis of stress–displacement relationships within specific cycles corresponding to 1st, 10th, 20th, and other fractional intervals of simulated results were conducted, as illustrated in Fig. 9. It is evident that alterations in the phase difference of tangential shear displacement significantly govern the shear stress–displacement response characteristics of the interface. Furthermore, shifts in both the initial shear locus and propagation direction (Fig. 2) were induced by modifications in the Y-direction displacement phase angle, resulting in distinct shear path morphologies. Pronounced discrepancies in the Y-direction stress response behavior (Fig. 9(b)) were consequently manifested through these kinematic transformations, which demonstrate the phase sensitivity of contact mechanics phenomena.

When the incremental phase shifts (δ_φ) = 0 or π (corresponding to D_p = π/2 and π/2 + π, respectively), the shear deformation trajectory at the contact interface exhibits an ideal circular pattern. Under this boundary condition, the unidirectional stability of shear direction during periodic loading ensures that the mechanical response curve maintains continuous nonlinear characteristics without abrupt discontinuities. When D_p = π/2 + π/2 and π/2 + 3π/2, the shear path of the interface degrades to a two-way linear pattern. Due to the reversal of shear direction during the shearing process, an abrupt change in the stress–displacement response relationship is induced when the shear displacement reaches the amplitude value.

With the exception of the initial cycle, the stress response curves in successive cycles exhibit substantial overlap, consistently displaying a conjugated stress path characterized by a bimodal response. This phenomenon becomes progressively more pronounced as the phase difference increases during the two transformation phases of D_p, specifically, when transitioning from π/2 to π/2 + π/2 and from π/2 + π to π/2 + 3π/2 (corresponding to a gradual evolution of the shear path from circular to bidirectional linear geometry). Conversely, during the subsequent transformations from π/2 + π/2 to π/2 + π and from π/2 + 3π/2 to π/2 + 2π (where the shear path transitions from bidirectional linear back to circular), the bimodal response diminishes in prominence as the phase difference grows.

4.2 Non-coaxial characteristics

Figure 10 illustrates the schematic correlation between the shear stress orientation and the shear displacement evolution during the 10th loading cycle, under conditions where D_p transitions from π/2 to π/2 + π/3 (refer to Experiments 2, 4–6). The shear mechanism reveals that the intrinsic tangential directionality of displacement propagation along the contact interface inherently aligns with the incremental shear displacement vector. However, the directional analysis demonstrates a persistent angular disparity between the shear stress trajectory and the local displacement increment direction. This systematic misalignment substantiates the existence of non-coaxial shear behavior within the contact plane during progressive shearing, which arises from the inherent phase-dependent mechanical response under coupled loading conditions.

The non-coaxial property constitutes a fundamental attribute of interface behavior under 3D cyclic shear loading scenarios, characterized by systematic quantification through non-coaxial angular responses. It is imperative to note that in 2D interface systems, the absence of stress–path coupling effects inherently precludes the manifestation of non-coaxial behavior. A comprehensive understanding of interfacial non-non-coaxial behavior enables precise determination of both the tangential deformation direction and its incremental evolution, thereby facilitating the formulation of a 3D numerical interface model with enhanced realism and predictive accuracy. This advanced modeling framework provides critical insights for investigating complex gravel–structure interaction mechanisms under multi-dimensional loading regimes.

The variation characteristics of the non-coaxial angle (α₃) at the interface were analyzed across 16 loading cases with distinct displacement phase differences. The results were presented in Fig. 11, with the following parameter definitions: α₁ represents the directional angle of tangential displacement increment, α₂ denotes the directional angle of tangential shear stress, and α₃ specifies the non-coaxial coupling angle between shear stress and displacement increments (α₃ = α₂ − α₁). Notably, when the phase difference D_p = π/2 + π/2 and π/2 + 3π/2, the shear deformation path transitions from a curvilinear profile to a bidirectional linear configuration. This phase-dependent geometric transition induces a partial loss of 3D shear characteristics at the interface, warranting exclusion from subsequent analysis. Statistical evaluation reveals that under all other phase difference conditions (excluding the aforementioned singular cases), α₃ exhibits pronounced periodic oscillations with enhanced magnitude. This persistent angular deviation quantitatively demonstrates substantial non-coaxial coupling effects between the shear stress orientation and displacement propagation direction during interface shearing processes.

When the displacement phase difference D_p = π/2 and π/2 + π, the interface exhibits a circular shear path. During the first loading cycle, the non-coaxial angle α₃ decreases monotonically from 90° to approximately 15° and stabilizes thereafter. This behavior underscores a synchronized relationship under circular shear paths: the 15° phase lag between α₁ and α₂ aligns with the classical shear lag theory, were displacement gradients lags behind shear stress propagation due to microstructural inertia and geometric constraints. Notably, as the shear path transitions from circular to linear (e.g., D_p increasing from π/2 to π/2 + π/2 or from π/2 + π to π/2 + 3π/2), the fluctuation amplitude of α₃ gradually amplifies. This augmentation is attributed to the asymmetric deformation of interfacial micro-convexities, where the geometric constraint of linear shear paths disrupts the equilibrium between microscale elastic interactions and macroscale kinematics, exacerbating directional misalignment between shear stress and displacement increment. Conversely, during the reversal from linear to circular shear paths (e.g., D_p decreasing from π/2 + π/2 to π/2 + π or from π/2 + 3π/2 to π/2 + 2π), the α₃ fluctuations attenuate, indicating a dominance of macroscopic slip mechanisms over micro-convexity-induced asymmetry. The restoration of circular geometry promotes symmetric stress redistribution, effectively mitigating directional deviations between shear stress and displacement increment. Compared to the fractional model of Sun et al. [31], the DEM simulations in this study capture the non-coaxial behavior induced by 3D shear paths, where the fluctuation amplitude of α₃ is directly related to the phase difference D_p.

It is underscored by these collective experimental findings that the regulatory control over the interfacial geometric constraint conditions is exerted by the displacement phase difference D_p, thereby dictating the coordination between the orientation of shear stress and the direction of the displacement increment. A path-dependent nature is demonstrated for the stress–strain relationship at the interface. The value of α₃ is shown not to be constant, but is predominantly governed by the geometry of the shear path, as quantified by D_p. The direction of plastic flow (increment of plastic strain) is indicated by the persistent non-coaxial behavior to be not solely determinable by the current stress state. A fundamental limitation of the conventional interfacial constitutive model is revealed to originate from the coaxiality assumption, which presumes the coincidence of the principal stress directions and the directions of the principal strain rates. The incorporation of an internal state variable representing the anisotropy of the microstructure is necessitated in future constitutive models for interfaces under 3D loading. The evolution law for this structural tensor is required to be coupled with the stress path, quantified by D_p, to enable accurate prediction of the non-coaxial angle α₃ and the associated energy dissipation under complex loading conditions.

4.3 Volumetric deformation

Figure 12 depicts the evolution of the normal displacement at the interface under shear conditions when D_p = π/2 + π/4. It is clearly evident that, as cyclic shearing progresses, the normal displacement at the interface demonstrates a gradual increase, which is accompanied by substantial volume compression of the gravel. However, over time, the rate of this increase progressively decelerates. The normal displacement (v) induced by cyclic shearing can be decomposed into irreversible and reversible components. These components correspond to the maximum magnitude (irreversible volumetric deformation is represented by black solid lines) and the recoverable component of the normal displacement within a single shear cycle, respectively. The irreversible normal displacement (v_ir) accumulates gradually, and this accumulation process slows down as the shear stress cycles. Throughout the evolution of the normal displacement, an obvious reversible normal volume change (v_re) emerged, as depicted in the shaded region of Fig. 12.

Figure 13 illustrates the evolution characteristics of the reversible normal volumetric deformation of the interface with tangential displacement under eight distinct phase difference conditions. Results from the 5th, 10th, 15th, and 20th cycles of each specimen’s cyclic shearing process were selected for detailed analysis. Observations reveal that as the shear displacement phase difference D_p evolves from π/2 to π/2 + π/2 (Figs. 13(a) and 13(e)) and transitions from π/2 + π to π/2 + 3π/2 (Figs. 13(c) and 13(g)) (where the shear path shifts from circular to bidirectional linear), the amplitude of the reversible normal displacement exhibits a gradual increasing trend. Conversely, when D_p progresses from π/2 + π/2 to π/2 + π (Figs. 13(b) and 13(f)) and from π/2 + 3π/2 to π/2 (Figs. 13(d) and 13(h)) (with the shear path transitioning from bidirectional linear back to circular), the amplitude of the reversible normal displacement demonstrates a progressive decreasing tendency. For cases with a circular shear path (D_p = π/2 and π/2 + π), the variation pattern of the reversible normal displacement remains less pronounced. This is attributed to the absence of abrupt changes in shear direction during the cyclic shearing process. In contrast, under bidirectional linear shear paths (D_p = π/2 + π/2 and π/2 + 3π/2), the reversible normal displacement develops significantly with progressive cyclic shearing, reaching its peak amplitude at the maximum shear displacement. This phenomenon arises because two shear direction reversals occur within each cycle during shearing. Notably, the peak amplitudes of the reversible normal displacement in positive and negative directions are nearly identical, indicating inherent isotropy in the interface’s volumetric deformation characteristics.

Figure 14 elucidates the evolution law of the reversible shear volumetric deformation on the interface as the displacement phase difference progresses from π/2 + π/6 to π/2 + π/2. Notably, when the shear path approximates a circular trajectory (with a phase difference D_p = π/2 + π/6), the peak value of the interface’s reversible normal displacement is markedly lower compared to the scenario where the shear path approaches bidirectional linearity (D_p = π/2 + π/3). This disparity is primarily attributed to the fact that controlling the phase difference of shear displacement results in a relatively minor abrupt change in the specimen’s shear direction, thereby impeding the full development of reversible normal displacement; conversely, a more significant abrupt change in the interfacial shear direction facilitates its complete development. Moreover, under displacement control conditions, with the continuation of cyclic shear, the peak value of the reversible normal displacement gradually diminishes before stabilizing. The maximum reversible normal displacement can reach approximately 2 mm and exhibits axial symmetry in both the x and y directions within a single shear cycle.

The observed phenomena reveals that the coupling between shear-path geometry and displacement phase difference significantly influences the thermodynamics of contact surface deformation. The insufficient development of reversible normal displacement under gentle shear-direction changes suggests limited energy dissipation through volume change, whereas abrupt directional changes enhance frictional interactions and associated volumetric responses. These findings provide a theoretical basis for predicting the deformation behavior of structural interfaces in geotechnical engineering, such as soil–structure interaction systems under complex cyclic loading conditions.

4.4 Microscopic responses

The contact force distribution at the beginning of the 10th cycle shear of the interface under different shear displacement phase differences when the normal stress σ_n = 400 kPa is shown in Fig. 15. Contacts are represented as solid lines connecting the centroids of two contacting particles, with line thicknesses proportional to contact force magnitudes. In the modeled system, soil particles were uniformly distributed within the shear box (Fig. 4(a)), which inherently leads to a relatively uniform horizontal contact force distribution within the specimen prior to cyclic shearing.

As illustrated, the red arrow denotes the initial shear direction of the specimen during its 10th loading cycle. Notably, variations of the tangential displacement phase difference not only alter the cyclic shear path of the interface but also modulate the initial shear orientation of the cyclic shearing process. Furthermore, with progressive cyclic shearing, the interparticle contact forces in the soil mass positioned rearward along the shear direction exhibit significantly greater magnitude and broader spatial distribution compared to those in the soil mass situated forward (distinguished by black dashed lines in Fig. 15). Consistent with the findings of Ari and Akbulut [28], the particle contact forces behind the shear direction are significantly enhanced. This observation indicates that the soil at the rear of the interface sustains greater shear resistance during shearing. This phenomenon arises from the application of shear forces to the trailing edge of the shear box’s movement direction, which promotes more effective interparticle force transmission within this rearward region.

Figure 16 presents the force chain distribution along the vertical cross-section of the specimen (oriented parallel to the shear direction) at the commencement of the 10th loading cycle under four distinct phase difference conditions. At the initiation of shearing, the contact force distribution among particles within the shear box is initially uniform. As shearing progresses, interparticle contacts intensify, leading to a dispersion of localized force concentrations along the specimen height. Notably, the lower region of the shear box exhibits a more homogeneous contact force distribution, with relatively smaller concentrated forces near the specimen base. In contrast, the upper portion of the shear box-closer to the structural panel-displays enhanced contact forces, forming a distinct inclined band extending from the rear to the front of the shear box. Within this band, interparticle contacts are significantly more pronounced, contrasting sharply with the lower shear box region. This discrepancy indicates the presence of a weak interaction zone with limited thickness between the gravel–structural units, approximately 5–7 times the median particle size D₅₀ of the soil (60 mm), as delineated in the red box in Fig. 16. Such a feature warrants careful consideration in both engineering design and safety assessment practices.

5 Conclusions

In this investigation, a 3D numerical model characterizing shear behavior at gravel–structure interface was developed based on DEM simulations. The influence of variations in the phase difference D_p of tangential shear displacement on the 3D cyclic shear response of the gravel–structure interface was examined. The manifestation of non-coaxial behavior during the shear process at gravel–structure interface were uncovered. Additionally, the microscopic mechanisms governing the effects of tangential shear displacement phase difference variations were elucidated through detailed investigations of the specimen’s internal structural evolution and shear kinematics. This study reveals a critical mechanistic insight: the geometry of the shear path, controlled by the phase difference D_p, is a fundamental governing parameter for the 3D cyclic response of gravel–structure interfaces, a factor previously unreported in the literature. The primary conclusions derived from this study are summarized as follows.

1) The critical role of shear path geometry controlled by the phase difference D_p was established in this study. It was demonstrated to be a fundamental governing factor for the interfacial mechanical response. The stress–strain hysteresis, the amplitude of non-coaxial fluctuations, and the development of reversible volumetric strain were all found to be critically determined by the shear path, moving beyond the conventional understanding based primarily on 2D loading.

2) Novel hysteretic and volumetric mechanisms were revealed. The hysteretic behavior was shown to be directly dictated by the shear path geometry. Stable, continuous hysteresis loops with minimal volumetric hysteresis were promoted by circular paths (D_p = π/2, 3π/2), whereas abrupt stress reversals and significantly amplified dilatancy-recompression cycles were induced by bidirectional linear paths (D_p = π, 2π). This provides a new mechanistic understanding of reversible volumetric deformation under complex 3D loading.

3) A path-dependent nature of non-coaxial behavior was discovered. Significant fluctuation of the non-coaxial angle (α₃), exceeding 30°, was observed under phase-shifted shear paths. The amplitude of these fluctuations was governed by phase difference D_p, with a transition from stable, damped oscillations under circular paths to large, periodic oscillations as the path evolved toward linearity. This discovery of path-dependent non-coaxial behavior is considered critical for advancing constitutive models beyond the co-axiality assumption.

4) A finite-thickness weak interaction zone at the interface was identified and quantitatively characterized from a micro-mechanical perspective. This zone measured approximately 5–7 times the mean particle diameter (D₅₀), was characterized by highly asymmetric force chain distributions and inclined stress propagation, thereby providing a physical basis for strain localization phenomena.

The explicit relationship between tangential displacement phase difference D_p, shear path geometry, and the resulting mechanical phenomena provides an essential benchmark for the development of advanced, path-dependent constitutive models for interfaces. Furthermore, the understanding of these mechanisms is deemed crucial for the improvement of stability analysis and seismic design of geotechnical structures subjected to multi-directional loading, ensuring a more accurate prediction of interface behavior in real-world conditions.

Despite the insights gained, several limitations of the proposed numerical approach are acknowledged. The mechanical behavior of irregular gravel was simulated using spherical particles with a rolling resistance contact model, which inherently simplifies the complex interlocking effects and fabric anisotropy dependent on actual particle morphology. Furthermore, the model’s validation under a specific set of conditions also warrants caution regarding its applicability to a broader range of gradations, densities, and complex loading paths. Addressing these limitations represents a promising direction for future research, aiming to enhance the generality and practical applicability of the numerical framework established in this work.

References

[1] Zhang C , Jiang C , Pang L . Investigation on the shear behavior of thin-layer grouted joints: Discrete element analysis and analytical model. Computational Particle Mechanics, 2025, 12(3): 1751–1771

[2] Yan Z , Zaoui A , Sekkal W . Discrete-element method study of the effect of ballast layer depth on the performance of railway ballast bed. International Journal of Geomechanics, 2025, 25(1): 04024309

[3] Afzali-Nejad A , Lashkari A , Martinez A . Stress–displacement response of sand–geosynthetic interfaces under different volume change boundary conditions. Journal of Geotechnical and Geoenvironmental Engineering, 2021, 147(8): 04021062

[4] Rehman Z U , Luo F , Wang T , Zhang G . Large-scale test study on the three-dimensional behavior of the gravel–concrete interface of a CFR dam. International Journal of Geomechanics, 2020, 20(6): 4020046

[5] Yang J , Li X S . State-dependent strength of sands from the perspective of unified modeling. Journal of Geotechnical and Geoenvironmental Engineering, 2004, 130(2): 186–198

[6] Jitsangiam P , Praai S , Boulon M , Jenck O , Chen X , Techavorasinsakul S . Characterization of a soil-rough structure interface using direct shear tests with varying cyclic amplitude and loading sequences under a large cyclic testing cycle condition. Acta Geotechnica, 2022, 17(5): 1829–1845

[7] Qannadizadeh A , Shourijeh P T , Lashkari A . Laboratory investigation and constitutive modeling of the mechanical behavior of sand–GRP interfaces. Acta Geotechnica, 2022, 17(10): 4253–4275

[8] Rui S , Wang L , Guo Z , Zhou W , Li Y . Cyclic behavior of interface shear between carbonate sand and steel. Acta Geotechnica, 2021, 16(1): 189–209

[9] Lashkari A , Jamali V . Global and local sand–geosynthetic interface behavior. Geotechnique, 2020, 71(4): 1–47

[10] Maghsoodi S , Cuisinier O , Masrouri F . Effect of temperature on the cyclic behavior of clay–structure interface. Journal of Geotechnical and Geoenvironmental Engineering, 2020, 146(10): 04020103

[11] Feng W , Li C , Yin J , Chen J , Liu K . Physical model study on the clay–sand interface without and with geotextile separator. Acta Geotechnica, 2019, 14(6): 2065–2081

[12] Sun T , Yang J , Shi K , Lin S , Yang J . Macro–microscopic mechanical study of clay–structure interface shear behavior using direct shear testing and DEM simulation. Advances in Civil Engineering, 2025, 2025(1): 6356879

[13] Feng D , Zhang J . Role of tangential displacement amplitude in 3D cyclic simple shear behavior of gravel–steel interface. International Journal of Geomechanics, 2024, 24(9): 04024189

[14] Feng D KZhang J M. Role of stress amplitude ratio in the non-coaxial behavior of the interface between gravel and structure. Rock and Soil Mechanics, 2022, 43(11): 3047–3058 (in Chinese)

[15] Feng D , Zhang J . Consistency behavior and mechanism of the irreversible dilatancy of gravel–structure interface subjected to 3D loadings. Acta Geotechnica, 2024, 19(7): 4637–4653

[16] Martinez A , Palumbo S , Todd B D . Bioinspiration for anisotropic load transfer at soil–structure interfaces. Journal of Geotechnical and Geoenvironmental Engineering, 2019, 145(10): 04019074

[17] Wang H , Zhou W , Yin Z , Jie X X . Effect of grain size distribution of sandy soil on shearing behaviors at soil–structure interface. Journal of Materials in Civil Engineering, 2019, 31(10): 04019238

[18] Wang J , Liu F Y , Wang P , Cai Y Q . Particle size effects on coarse soil-geogrid interface response in cyclic and post-cyclic direct shear tests. Geotextiles and Geomembranes, 2016, 44(6): 854–861

[19] Zhang G , Wang L , Zhang J . Dilatancy of the interface between a structure and gravelly soil. Geotechnique, 2011, 61(1): 75–84

[20] Feng D , Zhang J , Hou W . Three-dimensional stress-controlled cyclic behavior of the interface between gravel and steel. Soil Dynamics and Earthquake Engineering, 2021, 151: 106988

[21] Han F , Ganju E , Salgado R , Prezzi M . Effects of interface roughness, particle geometry, and gradation on the sand–steel interface friction angle. Journal of Geotechnical and Geoenvironmental Engineering, 2018, 144(12): 04018096

[22] Li M , Li Y , Islam M R . Effects of water content and interface roughness on the shear strength of silt–cement mortar interface. Soil and Foundation, 2021, 61(6): 1615–1629

[23] Chen J , Indraratna B , Vinod J S , Ngo T , Liu Y . Effects of particle shape on the shear behavior and breakage of ballast: A DEM approach. International Journal of Geomechanics, 2025, 25(1): 04024304

[24] Zeraati-Shamsabadi M , Sadrekarimi A . A DEM study on the effects of specimen and particle sizes on direct simple shear tests. Granular Matter, 2025, 27(2): 36

[25] Bernhardt M L , Biscontin G , O’Sullivan C . Experimental validation study of 3D direct simple shear DEM simulations. Soil and Foundation, 2016, 56(3): 336–347

[26] Zabihi NAthanasopoulos-Zekkos A. Discrete element modelling of large scale stacked-ring simple shear test of steel spheres. In: Proceedings of the Geo-Congress 2020: Modeling, Geomaterials, and Site Characterization. Minneapolis, MN: ASCE, 2020, 491–500

[27] Zhao J , Zhao S , Luding S . The role of particle shape in computational modelling of granular matter. Nature Reviews. Physics, 2023, 5(9): 505–525

[28] Ari A , Akbulut S . Evaluation of sand–geomembrane interface behavior using discrete element method. Granular Matter, 2022, 24(1): 21

[29] Stutz H H , Martinez A . Directionally dependent strength and dilatancy behavior of soil–structure interfaces. Acta Geotechnica, 2021, 16(9): 2805–2820

[30] DeJong J T , Westgate Z J . Role of initial state, material properties, and confinement condition on local and global soil–structure interface behavior. Journal of Geotechnical and Geoenvironmental Engineering, 2009, 135(11): 1646–1660

[31] Sun Y , Sumelka W , Gao Y , Nimbalkar S . Phenomenological fractional stress–dilatancy model for granular soil and soil–structure interface under monotonic and cyclic loads. Acta Geotechnica, 2021, 16(10): 3115–3132

[32] Sun Y , Wojciech S , He S , Gao Y . Enhanced fractional model for soil–structure interface considering 3D stress state and fabric effect. Journal of Engineering Mechanics, 2022, 148(9): 04022054

[33] Rehman Z U , Zhang G . Shear coupling effect of monotonic and cyclic behavior of the interface between steel and gravel. Canadian Geotechnical Journal, 2018, 56(6): 876–884

[34] Zia U R , Zhang G . Three-dimensional elasto-plastic damage model for gravelly soil–structure interface considering the shear coupling effect. Computers and Geotechnics, 2021, 129: 103868

[35] Rehman Z U , Zhang G . Cyclic behavior of the gravel-steel interface under varying rotational shear paths. Canadian Geotechnical Journal, 2019, 58(3): 305–316

[36] Liu Z , Gong G , Dou X , Zuo T , Zhao Y , Wu Q , Ning F . The microscopic failure mechanism of hydrate-bearing sand under direct shear: Insights from 3D DEM simulation. Energy & Fuels, 2025, 39(17): 8001–8010

[37] Liu Y , Liu H , Mao H . The influence of rolling resistance on the stress-dilatancy and fabric anisotropy of granular materials. Granular Matter, 2018, 20(1): 12

[38] Feng D , Zhang J , Hou W . Three-dimensional behaviors of a gravel–steel interface considering initial shear stress. Journal of Geotechnical and Geoenvironmental Engineering, 2021, 147(2): 04020158

RIGHTS & PERMISSIONS

Higher Education Press