Numerical Simulations of Extreme Deformation Problems in Granular-Dominated Hazard from Indoor to Engineering Geological Scale: A Comparative Study

1  Introduction

Granular material systems are widely present in both natural environments and industrial processes, and their dynamics often involve extreme deformation, multiphase coupling, and multiscale interactions. These flows exhibit inherent multiphase and multiphysics characteristics. Even a dry granular assembly can, depending on the local stress and shear rate, simultaneously exhibit different states that are analogous to solids, liquids, and gases. Furthermore, granular materials often interact with other interstitial fluids (e.g., water, air), giving rise to complex multi-phase flows such as debris flows (solid-fluid interaction) or gas-particle flows in industrial processes [1]. The dynamics of such systems involve strong couplings of momentum, mass, and energy across phases, posing significant challenges for numerical simulation. In industrial applications, granular flow problems range from micron-sized powder compaction and molding to particle transport in meter-scale hoppers. In geological hazards, large-scale phenomena such as landslides and avalanches occur at kilometer scales, frequently triggered by granular material instability. Therefore, investigating the multiscale large deformation mechanical behavior of granular systems holds significant scientific importance for understanding motion mechanisms, optimizing engineering equipment design, and supporting disaster prevention and mitigation.

Current numerical methods for granular flows fall into two primary theoretical frameworks: continuum methods and discrete methods [2]. Within the continuum approach, the motion of granular materials is governed by constitutive relations and described by solving continuum equations, implemented via grid-based techniques [3,4] or meshless particle methods. A wide range of constitutive relations are available, and the choice of a specific model depends on the flow regime and the physical mechanisms of interest [5].

Grid-based methods, such as the Finite Volume Method (FVM), Finite Element Method (FEM), and Finite Difference Method (FDM), discretize governing equations through integration or approximation over discrete cells or grid points. For geohazards where lateral dimensions dominate vertical scales, full 3D continuum solutions are often impractical. Consequently, depth-averaged models based on hydraulic principles have been developed [6–9]. These assume horizontal flow scale substantially exceeds vertical dimensions, satisfying shallow-water approximations, simplifying 3D equations to 2D systems closed with bed-stress models.

Meshless particle methods, such as Smoothed Particle Hydrodynamics (SPH) and the Material Point Method (MPM), approximate continuous fields through particle discretization. SPH, one of the earliest meshless techniques, originated in astrophysics and is now widely applied in the simulation of industrial granular flows [10] and small-scale landslides [11]. MPM combines the advantages of Eulerian and Lagrangian frameworks and is particularly effective for large-deformation problems in granular materials [12].

Discrete methods explicitly resolve particle contact mechanics, enabling precise dynamic analysis at the particle scale. While Molecular Dynamics (MD) effectively models nanoscale interactions, its computational cost limits large-scale applications. The Discrete Element Method (DEM), based on individual particle force analysis, studies granular assembly motion. DEM extensively simulates industrial processes and natural disasters, serving as a core tool for granular multiscale mechanics. The Discontinuous Deformation Analysis (DDA) method directly solves contact force-displacement relationships and is well-suited for applications such as rock fragmentation and pile-structure stability analysis.

In summary, the numerical simulation of large-deformation granular flows relies primarily on three categories of approaches: (1) discrete methods, which resolve individual particle interactions; (2) particle-based continuum methods, which discretize the continuum with Lagrangian particles; and (3) grid-based continuum methods. The last category includes efficient, simplified models such as the Depth-Averaged Model, which employs the shallow-flow assumption to reduce the full 3D equations to a 2D system, typically solved on a fixed grid (e.g., using FVM). Notably, the depth-averaging concept is versatile and can also be integrated with discretization methods [13]. Given the variety of available approaches, selecting the most appropriate method for a specific granular flow problem requires a clear understanding of their comparative strengths, limitations, and applicability across different scales.

To this end, this study systematically evaluates three representative methods, DAM, SPH, and DEM, by applying them to a suite of multiscale granular flow problems. These range from granular column collapse at the particle scale, to flow-structure interaction at the laboratory scale, and the 2015 Shenzhen Guangming landslide at the field scale. While each method offers strong capabilities for simulating extreme deformations, differences in their physical assumptions, governing equations, and discretization approaches fundamentally affect their predictive accuracy, computational efficiency, and thus their engineering suitability. Our evaluation involves a multi-level verification test suite to assess operational conditions, computational performance, and scale applicability.

2  Numerical Model

2.1 Depth-Averaged Model (DAM)

Savage and Hutter [6,7] established the depth-averaged model for thin-layer granular flows by combining the shallow-flow assumption with Coulomb’s friction law, resulting in a set of hyperbolic governing equations. This equation set exhibits mathematical similarities to classical shallow-water equations. The resulting model, termed the Savage-Hutter (SH) equations, has been extensively studied and extended [8,14,15].

Regarding numerical solutions for the SH model, early studies employed Lagrangian finite-difference methods constrained by grid distortion. With advances in computational fluid dynamics, Eulerian Godunov-type methods became mainstream [16]. To capture shock waves, advanced algorithms like oscillation-free schemes and central difference schemes were applied in the early 21st century [17]. Pudasaini et al. [18] implemented an inherent oscillation-free scheme for the SH model, validated through laboratory-scale thin-layer granular flow simulations. These improvements enhanced computational accuracy and facilitated widespread engineering applications in landslide and debris flow simulations.

2.1.1 Convergence Analysis

In numerical computation, grid-based methods exhibit rigorous mathematical convergence. Convergence describes how discrete solutions approach continuous solutions under grid refinement or increased iterations. Two essential conditions are consistency, meaning the discrete equations approximate the differential equations, and numerical stability, which requires controlled error propagation. The Lax equivalence theorem states that for well-posed linear problems, consistency and stability guarantee convergence [19].

For example, solving the single-wave equation via second-order finite volume methods (FVM) achieves a convergence order of 1.9 in smooth regions, which closely matches the theoretical second-order accuracy [20]. Shock waves or discontinuities induce numerical oscillations due to conservation law nonlinearity. Artificial viscosity stabilizes computations but reduces accuracy to first-order [17]. This trade-off is fundamental: Godunov’s theorem states that no linear scheme maintains monotonicity and higher-than-first-order accuracy at discontinuities [21]. High-order weighted essentially non-oscillatory (WENO) schemes address this by adaptively adjusting stencil weights, thus preserving high-order accuracy in smooth regions while reducing order near discontinuities [22,23].

2.1.2 Boundary Conditions

Boundary condition treatment is a core advantage of grid-based methods over meshless approaches. The Finite Volume Method (FVM) and Finite Element Method (FEM) rigorously enforce Dirichlet (Type I), Neumann (Type II), and Robin (Type III) boundary conditions [21]. For instance, in fluid dynamics, no-slip conditions at solid boundaries are implemented through Dirichlet conditions that specify zero velocity, while convective heat transfer relies on Robin conditions involving mixed derivatives. This direct physical-mathematical correspondence provides a robust foundation for practical applications.

2.1.3 Physical Model and Numerical Method Limitations

Continuum methods require constitutive models to establish macroscopic stress–strain relationships and close governing equations. Granular flow constitutive models include: phenomenological models based on steady-state flow characteristics [5], notably the widely applied µ(I)-rheology model [24–26], and kinetic theory incorporating particle size effects and inelastic collisions [27,28]. Despite diversity, model selection and parameter calibration remain challenging. A key continuum method limitation thus lies in constitutive model validity. Depth-averaged models partially mitigate this by simplifying stress models via the shallow-flow assumption [29].

Numerically, FEM and FVM differ in formulation: FEM solves algebraic systems from differential equation weak forms, while FVM uses conservation law integral forms. FEM discretizes domains into elements, approximates solutions via shape functions, and employs the Galerkin weighted residual method to derive stiffness or mass matrix systems. However, weak formulations lack guaranteed local conservation, risking mass/momentum accumulation, while storing and solving large stiffness matrices demands significant memory. Conversely, FVM partitions domains into non-overlapping control volumes, integrates conservation equations per volume, and applies Gauss’s divergence theorem to convert volume integrals to surface fluxes. Core steps involve flux computation, time discretization, and source term linearization [30].

2.2 Smooth Particle Hydrodynamics (SPH)

The Smooth Particle Hydrodynamics (SPH) method originated in 1977 when Lucy and Monaghan proposed it for non-axisymmetric astrophysical gas dynamics [31,32]. This meshless Lagrangian approach discretizes continuous media into particles that carry physical properties, such as mass, density, and velocity. Field functions and derivatives are approximated through kernel-based weighted interpolation, while particle evolution is governed by conservation equations.

When extended to granular systems, SPH treats granular flows as non-Newtonian fluids, closing governing equations via stress models. Adaptive smoothing length techniques improve computational efficiency by dynamically adjusting particle influence domains [33]. GPU-optimized parallel neighbor search algorithms accelerated large-scale particle computations [34].

2.2.1 Convergence Analysis

SPH convergence lacks a rigorous mathematical establishment, with limited related studies. In 2003, Li’s team analyzed particle method convergence using two weighting functions, demonstrating sub-first-order accuracy [35]. Oger et al. [36] and Xu et al. [37] showed that uniform particle distribution with ideal kernels achieves second-order accuracy in interior regions but only zero-order accuracy at boundaries. Non-uniform particle distributions reduce interior accuracy below first-order. Despite lacking unified theory, case-specific particle resolution verification persists [34].

Higher-order consistency corrections address accuracy limitations, including: Conservative Reproducing Kernel SPH (CRKSPH) enabling linear field interpolation while conserving mass, momentum, and energy; Renormalized SPH (RSPH) and Modified SPH (MSPH) restoring kernel consistency via correction matrices; and Incompressible SPH (ISPH) eliminating oscillatory errors through pressure Poisson solutions [38].

2.2.2 Boundary Conditions

SPH handles boundary conditions with greater complexity than grid-based methods. To prevent particle penetration at solid walls, techniques such as the mirror particle method [39] and the repulsive force method [40] are employed. The former generates virtual particles for kernel interpolation, while the latter applies empirical potentials near the walls. These methods enforce free-slip or no-slip conditions, corresponding to Dirichlet (Type I) boundaries. However, SPH lacks direct support for Neumann (Type II) or Robin (Type III) conditions. In Dissipation-Free SPH (DFSPH), pressure gradients are implicitly imposed by solving the constraint ∇⋅v=0, which significantly increases computational cost [41].

2.2.3 Physical Model and Numerical Method Limitations

As a continuum approach, SPH requires constitutive models to close stress equations. Particle-based discretization introduces challenges, including irregular particle distributions degrading kernel approximations and inducing strain-rate errors, and boundary particle deficiency causing inaccurate wall constitutive responses [42].

SPH typically employs explicit integration schemes: second-order leapfrog (LF), predictor-corrector, and Runge-Kutta (RK) methods. The LF algorithm is favored for low memory requirements and single-step optimization per iteration.

2.3 Discrete Element Method (DEM)

The Discrete Element Method (DEM) was introduced by Cundall and Strack [43]. This method captures microscopic mechanical behavior in granular systems by tracking individual particle motion and interactions. Its core concept treats particles as independent entities, modeling instantaneous collisions and friction through contact mechanics. DEM circumvents traditional continuum theories’ reliance on macroscopic constitutive relationships, directly capturing particle trajectories. It excels in simulating granular accumulation, flow evolution, and shear-induced failure [42,44].

However, real-world geotechnical problems involve millions to billions of particles. Direct microscopic interaction simulation incurs exponentially increasing computational costs. Coarse-graining techniques alleviate this bottleneck by grouping clusters of microscopic particles into equivalent macroscopic super-particles while preserving original statistical mechanical properties through careful parameter calibration [45,46]. Although this approach enhances computational efficiency, the success of coarse-graining critically depends on accurate micro-to-macro parameter mapping. This process determines the equivalent mechanical properties of super-particles through DEM simulations or experimental data, but inaccurate calibration may distort macroscopic flow characteristics [47].

2.3.1 Convergence Analysis

DEM lacks mathematically rigorous convergence concepts typical of continuum mechanics. Instead, convergence manifests through macroscopic statistical behavior approximation via accurate microscopic interaction modeling, rather than strict mathematical criteria.

At micro-to-meso scales, DEM directly assigns physically meaningful parameters consistent with real particles, avoiding theoretical convergence discussions. At larger scales, computational constraints necessitate coarse-graining, replacing real particles with size-distributed super-particles. Here, DEM convergence refers to stabilization of system-level statistical quantities such as stress tensors under increasing simulated particle counts.

2.3.2 Boundary Conditions

DEM simulates geometric constraints using wall particles or imported entities, effectively implementing Dirichlet-type displacement boundaries. Without continuum field-variable constructs, DEM cannot directly implement Neumann- or Robin-type conditions.

2.3.3 Physical Model and Numerical Method Limitations

DEM fundamentally differs from continuum mechanics by constructing particle system dynamics directly from microscopic contact mechanics, not constitutive stress–strain relationships. Contact model refinements extend DEM to complex granular systems [42].

DEM employs explicit integration for coupled ordinary differential equations governing particle motion. Each time step initiates with contact detection. Traditional brute-force methods exhibit O(N2) time complexity, while spatial partitioning or hierarchical data structures achieve near-linear complexity [48]. After contact force computation, motion equations are solved via explicit integration: single-step [43], multi-step [49], and predictor-corrector methods. Integration consumes minimal time per step, but algorithm selection critically determines maximum allowable time steps for large-scale simulation acceleration.

2.4 Comparative Summary of the Three Methods

The comparative analysis of DEM, SPH, and DAM reveals fundamental distinctions arising from their theoretical formulations and discretization frameworks. DEM is rooted in a discrete description, explicitly resolving interparticle contacts by solving Newton’s laws of motion for each particle. In contrast, SPH and DAM are both continuum-based but adopt different discretization approaches: SPH utilizes moving particles to solve the full three-dimensional governing equations, while DAM first applies the shallow-flow assumption to reduce the physics to a two-dimensional system, which is then typically discretized on a fixed grid. These foundational differences define their respective capabilities, limitations, and natural domains of applicability across scales, as systematically summarized in Table 1.

The comparison underscores a fundamental trade-off between physical fidelity and computational efficiency, dictated by the chosen modeling framework. DEM provides the most direct physical representation at the grain scale but is constrained by computational cost. SPH offers a balance, capturing complex deformations in a continuum framework but requiring careful handling of constitutive models and boundaries. DAM achieves high efficiency for large-scale problems by leveraging a simplifying physical assumption (shallow flow), which in turn defines its applicability boundary. Thus, the selection of an appropriate method is not merely a technical decision but is fundamentally guided by the scale of interest, the required physical detail, and the available computational resources.

3  Multiscale Case Study Setup

Validation tests represent specialized computational fluid dynamics (CFD) experiments designed for model verification, requiring precise accuracy specifications. Employing a hierarchical validation framework to verify complex systems, these tests constitute the most reliable model confirmation approach [50]. This framework decomposes system-level validation into four progressive tiers through a layered strategy: unit problem, benchmark problem, subsystem problem, and complete system problem (Fig. 1a). Geometric and flow condition complexity escalates across tiers, following a validation path from simple isolated systems to complex coupled systems.

Figure 1: (a) Hierarchical validation structure for confirmation tests. (b) Verification levels corresponding to characteristic granular flow scales.

To systematically evaluate the reliability of the three numerical methods for simulating multiscale and multiphysics granular flows, this study establishes a layered validation framework based on confirmation test theory. Through comparative analyses using DEM, SPH, and DAM on selected benchmark cases, the research delineates each method’s applicability boundaries in capturing macroscopic flow behavior. The numerical implementations are as follows. DEM simulations were conducted using the open-source LIGGGHTS platform, utilizing a soft-sphere approach with the Hertz–Mindlin contact model [51]. SPH simulations were based on the DualSPHysics code, extended with in-house developments for granular materials. An elastoplastic constitutive model based on the Drucker–Prager yield criterion was employed, stabilized by artificial viscosity and a repulsive-force boundary condition [34,52]; the smoothing length was set as h = dp, with dp being the initial particle spacing. DAM simulations were performed using a high-performance in-house code based on a finite-volume framework [8,9]. It implemented the Savage–Hutter depth-averaged model grounded in the Mohr–Coulomb criterion, utilizing a high-resolution Kurganov-type Riemann solver to capture flow fronts and discontinuities [53]. Case-specific parameters, including particle counts (for DEM/SPH) and grid resolutions (for DAM), are detailed in the corresponding setups (Table 2).

The unit problem addresses fundamental physical law representation fidelity, for example, energy conservation in particle contact mechanics, extensively validated by classical studies [43,44]. Validation herein therefore commences at the benchmark level.

Selected benchmark problems include classical 2D and 3D dam-break cases. Using controllable initial conditions and standardized geometries, these problems validate the method's accuracy in capturing mesoscale flow physics. The subsystem problem introduces granular flow-structure interactions to evaluate method performance under complex laboratory-scale boundaries. The complete system problem constructs a kilometer-scale geological model based on the 2015 Shenzhen landslide, assessing each method’s terrain adaptability and computational efficiency in large-scale geological flows to inform engineering disaster prediction method selection.

3.1 Benchmark Problem: Granular Dam Break

The granular dam-break problem serves as a benchmark case in granular flow simulations due to its reproducible geometry, representative flow mechanisms, and accessible high-quality experimental data. This classical setup features standardized geometric dimensions, including initial pile height and width, and controlled release conditions such as instantaneous gate opening, providing a repeatable physical scenario for numerical validation.

3.1.1 Case Setup

Simulation data reference laboratory experiments by Lajeunesse et al. [54], who investigated glass bead collapse in rectangular and cylindrical channels on horizontal planes. The study analyzes granular pile flow behavior with varying height-to-width ratios.

The quasi-2D setup features a rectangular channel (100cm×30cm×45mm) with a vertical gate (Fig. 2a). The initial pile geometry is defined by the gate position and initial height. The narrow channel width ensures unidirectional flow upon gate release. The cylindrical granular pile collapse represents a three-dimensional axisymmetric flow process (Fig. 2b). The configuration includes a cylindrical retaining gate and horizontal base, with initial pile geometry defined by height. Experimental constraints necessitated a half-cylinder gate in the original experiment; however, the numerical simulation adopts a full-cylinder model to eliminate boundary effects and ensure axisymmetry. Initial height-to-width ratios follow Lajeunesse et al.’s experimental configurations. Three aspect ratios represent distinct flow regimes.

Figure 2: (a) Quasi-2D collapse experimental setups. (b) 3D axisymmetric collapse experimental setups.

The particle parameters in experiments and numerical simulations are summarized in Table 2.

3.1.2 2D Granular Dam Break Results

Extensive research on granular dam-break flow characteristics with varying height-to-width ratios demonstrates that initial geometric parameters critically influence flow patterns and pile morphology. In 2D granular dam-break processes, lateral expansion along the rectangular channel transitions final pile shapes from trapezoidal to triangular (Fig. 3e). This phenomenon closely relates to internal granular flow dynamics. This section compares three numerical methods’ simulation results under different height-to-width ratios, analyzing consistency with experimental observations and method limitations.

Figure 3: Dynamic changes in the particle pile profile during the 2D granular dam-break process.

For a low height-to-width ratio (a = 0.6), the granular pile exhibits a trapezoidal shape upon flow cessation. Post-gate opening, the system displays layered motion: bottom particles move first, forming a rapid shear layer, while top portions near sidewalls remain quasi-static until flow stops. Numerical simulations show all three methods predict pile height and motion distance, agreeing well with experiments. DEM and SPH accurately capture flow dynamics and front progression details. However, DAM satisfies shallow-water conditions, causing bottom-first motion initiation and inaccurate transient particle morphology predictions.

At an increased height-to-width ratio (a = 2.4), flow patterns change. Post-release, the pile undergoes subsidence with a flat top surface, followed by bottom channel extension and complete top collapse, forming a triangular pile. DEM and SPH maintain high accuracy, but DAM underpredicts the final sediment height. Non-negligible vertical velocity/acceleration generates non-hydrostatic pressures within the granular column, violating shallow-water assumptions and causing DAM deviations.

Under a high height-to-width ratio (a = 16.7), strong non-hydrostatic characteristics emerge. The granular column experiences initial free-fall, tilting toward the channel, and final triangular pile formation. DEM and SPH predict sediment height and motion distance matching experimental data. However, the experimental central depression contrasts sharply with simulated convex profiles. SPH exhibits non-physical top particle-sidewall separation during settlement due to boundary flaws. Extreme conditions invalidate hydrostatic assumptions, rendering DAM ineffective.

3.1.3 3D Granular Dam Break Results

The 3D granular dam-break exhibits a flow regime transition analogous to the 2D case, governed by the initial height-to-width ratio. All three methods effectively captured the macroscopic dynamics and final runout across the tested ratios (Fig. 4). However, a key distinction from the 2D simulation was the pronounced lateral expansion in 3D, which enhanced inter-particle contacts and frictional dissipation. This led to greater deviations in predicting the final deposit height compared to 2D for all methods.

Figure 4: Dynamic changes of cross-sectional pile profiles during the 3D granular dam-break process.

The comparative performance of the methods in 3D reinforced the trends observed in 2D, but with amplified computational implications. DAM remained efficient for low aspect ratios but was similarly limited by its shallow-flow assumptions. SPH handled the large deformation and free surface well, though its boundary treatment and constitutive model limitations persisted. DEM continued to deliver the most accurate results by explicitly resolving 3D contact mechanics and energy dissipation, but the computational cost escalated dramatically due to the increased number of particles required for a 3D simulation, underscoring its primary bottleneck for large-scale problems.

3.2 Subsystem Problem: Particle Flow and Structure Interaction

Following benchmark verification under simple geometric conditions, this study advances to a more complex subsystem-level validation. This stage addresses engineering scenarios with complex boundaries, including granular material interactions with baffles during hopper discharge, debris flow impacts on bridge piers, and avalanche shocks on protective structures. These simulations involve highly complex granular flow states where particle-structure interactions intensify problem difficulty. Laboratory-scale particle flow-structure interaction at meter-to-decameter scales realistically reproduces such phenomena, providing an ideal platform for validating numerical method applicability [55].

3.2.1 Example Setup and Elevation Change at Measurement Points

This section simulates a typical particle flow-structure interaction experiment by Juez et al. [56], wherein granular piles flow along inclined planes interacting with three downstream hemispherical obstacles. Fig. 5a shows the model schematic: the left shaded area indicates the initial hemispherical granular pile, while the right displays hemispherical obstacles. Four measurement points around obstacles extract granular flow depth for tri-method comparison.

Figure 5: (a) Terrain setup for granular pile impact on hemispherical structures on a slope. (b) Elevation change process at measurement points.

Absent specific bed friction angles in literature, rigorous parameter calibration established that setting both bed friction angle and particle internal friction angle to 34° enables accurate reproduction of experimental pile morphology and spreading characteristics [53]. Simulation physical parameters are listed in Table 2.

Four measurement points in geometrically sensitive obstacle regions, specifically at the large ball shoulder and sides, compare experimental and simulated depth-time evolution curves, analyzing method performance in complex flows.

Fig. 5b shows that for points PS0L and PS0R, DAM most accurately predicts the depth-change trends. DEM captures an experimental peak when the flow reaches the points and later converges to stable depths consistent with other methods. SPH shows significant deviation in estimating flow arrival time, but eventually converges to stable depths similar to those of DEM and DAM. For PSL or PSR, DEM and DAM predict similar peak depths and stable depths matching experiments. SPH significantly overpredicts flow depth, indicating reduced reliability in complex boundary precision.

The deviations produced by the three methods are mainly attributed to the following factors: DAM’s shallow-water assumptions conflicting with significant 3D effects near obstacles, reducing momentum transfer, and weakening shock waves at side points while overpredicting shoulder-point peaks. DEM accurately simulates interactions, but yields elevated statistical values due to particle discretization. SPH’s boundary limitations cause large prediction errors in complex boundaries, compromising result precision.

3.2.2 Flow State Comparison

Granular flow-obstacle interaction involves complex dynamics. Upon impact, particles bypass obstacles, substantially altering flow depth, velocity, and direction. As shown in Fig. 6, at t = 460 ms, the flow front approaches the smaller obstacle. Contact with the larger hemisphere occurs at t = 500 ms. Direct impact at t = 640 ms generates three intensity-varied shock waves. These gradually merge into a single wavy shock by t = 740 ms, which relaxes after t = 900 ms, ultimately forming stable deposits. Flow stabilizes completely by t = 1500 ms.

Figure 6: Granular flow impact on three hemispherical structures during flow state evolution (a) Experiment, (b) DAM, (c) SPH, (d) DEM.

Fig. 6 reveals method-dependent flow evolution and deposition differences. DEM excels in reproducing dynamics, accurately capturing bypassing, climbing, and shock wave evolution. Despite edge particle dispersion, overall flow characteristics align closely with experiments. SPH results exhibit excessive flow aggregation, failing to reproduce experimental tail diffusion and oversimplifying obstacle wake structures. This stems from spatial averaging of interparticle interactions, weakening discrete collision characteristics and reducing complex-boundary flow resolution. DAM effectively predicts overall flow processes and final deposition but deviates in bypassing motion and vertical climbing behavior during particle-structure interaction. Given minimal deposition thickness in wake regions, results remain broadly consistent with experiments.

All three methods show high consistency in overall evolution, capturing climbing or bypassing phenomena, obstacle-front shock waves, and rear vacuum regions. However, each exhibits varying deviations in flow detail resolution.

3.3 Complete System Problem: Shenzhen Landslide Disaster

The final-tier system problem addresses engineering geological scales. Using Shenzhen’s 2015 Guangming New District landslide as an engineering case, a kilometer-scale geological model simulates large-scale disaster prediction, comparing method performance. High-speed long-runout landslides feature large volumes, long travel distances, and strong terrain adaptability [57]. Their dynamics involve multi-field coupling, complex topography, and cross-scale evolution, presenting dual challenges to numerical method precision and efficiency. This section evaluates DEM, SPH, and DAM engineering applicability regarding granular flow propagation along complex terrain and deposition pattern prediction.

Surface and initial source topography elevation data are derived from Ouyang et al. [3]. DAM parameters follow literature, implementing the Mohr-Manning model with parameter verification, where the pore pressure coefficient is 0.75, and the Manning coefficient is 0.01. SPH parameters reference previous work [34,52], where cohesion is 500 (Pa). DEM parameters align with SPH settings, and the specific parameter settings are shown in Table 2.

Simulation results (Fig. 7) show the red line indicating Ouyang et al.’s disaster boundary. All methods’ mainstream flow domains remain within the measured boundary, with flow front trajectories broadly matching field observations, particularly near building interactions. DEM simulations exhibit particles extending beyond the boundary, likely due to DEM parameter settings assigning equivalent particle sizes over tenfold larger than actual conditions. This increases flow discretization, causing deviations from overall trends. Alternatively, computational terrain exclusion of the original destroyed buildings in flow zones may contribute to discrepancies.

Figure 7: Flow state changes in the Shenzhen landslide (a) DAM, (b) SPH, (c) DEM.

Quantitative final disaster boundary analysis (Table 3) shows that all methods predict disaster length and width within 20% error of actual values. All three methods thus demonstrate applicability for large-scale geological disaster prediction precision. However, method selection warrants consideration beyond apparent accuracy due to parameter calibration difficulty and computational efficiency differences.

4  Analysis and Discussion

This section systematically analyzes DEM, SPH, and DAM adaptability across scales, considering usage conditions, applicable ranges, and computational efficiency.

4.1 Applicable Range and Usage Conditions

In granular material research, the four verification levels of the confirmation test framework correspond to characteristic granular flow scales (Fig. 1b).

Unit Problem (Particle Scale): Continuum methods (DAM, SPH) cannot resolve particle-scale contact behavior due to theoretical limitations. DEM explicitly solves Newton’s equations of motion by directly assigning physical parameters such as elastic modulus and friction angle, which makes it well-suited for simulation at this scale.

Benchmark Problem (Representative Elementary Volume, REV Scale): DAM effectively reproduces macroscopic shapes such as those in dam-break scenarios, yet it is limited by shallow-flow assumptions and cannot capture mesoscopic parameters like shear rate or concentration. SPH captures mesoscopic flow fields but is prone to stress instability and relies heavily on constitutive models. DEM accurately simulates granular dynamics and establishes stress–strain relationships, maintaining its effectiveness despite being sensitive to rolling friction parameters.

Subsystem Problem (Laboratory Scale): DAM efficiently predicts flow states with convenient calibration, but it simplifies vertical motion and lacks the ability to calculate structural loads. SPH effectively captures free surfaces, yet it exhibits boundary anomalies such as sidewall particle separation and struggles with multiphase transitions. DEM captures complex flow states during interactions, though the use of simplified particle shapes necessitates error mitigation through calibration.

Complete System Problem (Engineering Geological Scale): DAM predicts kilometer-scale flow ranges, deposition patterns, and velocity evolution with minimal parameters, enabling rapid disaster prediction. SPH models complex terrain and large deformations, but requires constitutive parameter validation when scaling from meters to kilometers. DEM faces prohibitive computational costs in simulating realistic particle distributions. Although coarse-graining by enlarging equivalent particles enables large-scale simulation, it introduces convergence issues and necessitates careful verification of scaled parameter correspondence.

Table 4 summarizes the three methods’ applicability across scales based on this analysis.

4.2 Computational Efficiency

Beyond applicability and accuracy, computational efficiency constitutes a crucial consideration when selecting numerical methods for multiscale granular flow problems. This section compares actual computation times for DAM, SPH, and DEM across the three scale cases simulated in Chapter 3. We acknowledge that a fully controlled hardware benchmark is ideal for such comparisons [58]. While constrained by available computational resources, we have standardized key parallel configurations to mitigate variability: DEM simulations were executed on an Intel Core i7-11700K processor, while SPH and DAM were run on an Intel Core i7-10700 processor. Critically, all three methods were configured to utilize an identical setup of 4 physical cores and 8 threads during parallel computation. This approach allows us to focus on revealing the inherent algorithmic scaling trends.

As detailed in Table 5, DAM consistently demonstrates the shortest computation times across all problem scales. SPH exhibits intermediate efficiency, while DEM requires the longest durations for dam-break and granular flow impact simulations, though it remains within acceptable limits. For large-scale landslides, however, DEM incurs prohibitively high computation times relative to other methods. These observed trends are a direct consequence of the fundamental algorithmic complexities: DEM’s explicit contact resolution scales with particle number, SPH’s particle-based continuum solver incurs neighbor search costs, and DAM’s depth-averaged equations are solved efficiently on a fixed grid. Without simplifications like coarse-graining, DEM’s computational costs become prohibitive for engineering applications at field scales.

5  Conclusion

This study verifies the applicability of three typical numerical methods for extreme deformation and flow problems in granular material systems. Three representative problems, including dam-break (element scale), granular flow impact on structures (laboratory scale), and the Shenzhen landslide (engineering scale), are simulated by DEM, SPH, and DAM methods for comparisons.

All three methods are capable of capturing key macroscopic dynamics, such as flow propagation and final deposition morphology, and provide quantitative metrics like runout distance. However, their performance diverges significantly due to fundamental differences in their underlying principles. At the meso- to macro-scales, DEM, by explicitly resolving particle-scale contact mechanics, accurately reproduces shear banding and energy dissipation, albeit at a computational cost that scales with particle number. At the laboratory scale, SPH shows strong adaptability to large deformations and free-surface evolution but is challenged by boundary treatment issues and simplified constitutive models, which can lead to non-physical particle separation and errors in local flow depth. For engineering-scale problems, DAM provides an efficient prediction of flow extent and deposition patterns using minimal parameters, yet its reliance on the shallow-flow assumption limits its accuracy in flows with strong non-hydrostatic effects or significant vertical momentum, and it cannot directly provide structural load estimates.

Regarding method selection, DEM is recommended for investigations focusing on fundamental granular mechanics or high-fidelity analysis of complex flow-structure interactions, especially when accelerated by parallel computing. SPH represents a balanced choice for problems prioritizing computational efficiency while retaining the ability to model large deformations and free surfaces, provided that boundary algorithms are improved to mitigate non-physical separation. DAM is the most suitable tool for rapid hazard assessment and large-scale geological flow simulations where efficiency is paramount; incorporating terrain gradient corrections within non-hydrostatic pressure terms can improve its accuracy on steep slopes.

For integrated cross-scale analysis, a hybrid strategy is promising. Domain decomposition allows coupling DEM (for near-field, mechanics-critical zones) with SPH or DAM (for far-field, large-deformation regions) to balance accuracy and efficiency. Alternatively, parameter transfer frameworks can employ microscopic DEM simulations to calibrate the rheological parameters of continuum models like DAM, thereby establishing a micro-to-macro validation chain that overcomes the scale limitations of individual methods.

Beyond the methodological insights discussed above, this study also has several limitations that point to directions for future research. First, the physical complexity was intentionally simplified; critical effects in many natural granular flows, such as water-particle coupling, particle size segregation, and basal erosion, were not considered. Second, the comparative evaluation was confined to three representative methods within distinct frameworks. Future work could broaden this scope. Incorporating other advanced approaches like the Particle Finite Element Method (PFEM) or the Material Point Method (MPM) would map a more comprehensive landscape of available simulation tools. Finally, while efforts were made to ensure a fair comparison of computational efficiency, a strictly controlled benchmark using unified hardware and software platforms remains an ideal goal for future method evaluations. Addressing these aspects will further advance robust, multi-scale numerical simulation capabilities for granular-dominated geohazards.

Acknowledgement: Not applicable.

Funding Statement: This work is financially supported by the National Natural Science Foundation of China (Nos. 12572465, 12032005).

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Yuxin Tian, Xiaoliang Wang; data collection: Yuxin Tian, Wanqing Yuan, Wangxin Yu; analysis and interpretation of results: Yuxin Tian, Xiaoliang Wang; draft manuscript preparation: Yuxin Tian; management and coordination: Xiaoliang Wang, Qingquan Liu; funding: Xiaoliang Wang, Qingquan Liu. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that support the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: None.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Ghannam A, Chehade A, Generous MM, Alazzam A, Kleinstreuer C, Ahmadi G, et al. A comprehensive review of particle-laden flows modeling: single/multiphase modeling approaches, benchmarks, heat transfer, intermolecular interactions, recent advances and future directions. Phys Rep. 2025;1118(4):1–96. doi:10.1016/j.physrep.2025.03.001. [Google Scholar] [CrossRef]

2. Hogue C, Newland D. Efficient computer simulation of moving granular particles. Powder Technol. 1994;78(1):51–66. doi:10.1016/0032-5910(93)02748-Y. [Google Scholar] [CrossRef]

3. Ouyang C, Zhou K, Xu Q, Yin J, Peng D, Wang D, et al. Dynamic analysis and numerical modeling of the 2015 catastrophic landslide of the construction waste landfill at Guangming, Shenzhen, China. Landslides. 2017;14(2):705–18. doi:10.1007/s10346-016-0764-9. [Google Scholar] [CrossRef]

4. Li G, Lin X, Lv Z, Wang Y, Zeng K, Wang F, et al. Continuous particle method for granular particle flow simulation and its extended application. AIP Adv. 2025;15(8):085223. doi:10.1063/5.0231069. [Google Scholar] [CrossRef]

5. Forterre Y, Pouliquen O. Flows of dense granular media. Annu Rev Fluid Mech. 2008;40(1):1–24. doi:10.1146/annurev.fluid.40.111406.102142. [Google Scholar] [CrossRef]

6. Savage SB, Hutter K. The motion of a finite mass of granular material down a rough incline. J Fluid Mech. 1989;199:177–215. doi:10.1017/s0022112089000340. [Google Scholar] [CrossRef]

7. Savage SB, Hutter K. The dynamics of avalanches of granular materials from initiation to runout. Part I: analysis. Acta Mech. 1991;86(1):201–23. doi:10.1007/BF01175958. [Google Scholar] [CrossRef]

8. Wang X, Li J. A new solver for granular avalanche simulation: indoor experiment verification and field scale case study. Sci China Phys Mech Astron. 2017;60(12):124712. doi:10.1007/s11433-017-9093-y. [Google Scholar] [CrossRef]

9. Wang X, Liu Q. Modeling shallow geological flows on steep terrains using a specific differential transformation. Acta Mech. 2021;232(6):2379–94. doi:10.1007/s00707-021-02944-3. [Google Scholar] [CrossRef]

10. Bui HH, Nguyen GD. Smoothed particle hydrodynamics (SPH) and its applications in geomechanics: from solid fracture to granular behaviour and multiphase flows in porous media. Comput Geotech. 2021;138(5):104315. doi:10.1016/j.compgeo.2021.104315. [Google Scholar] [CrossRef]

11. Ikari H, Gotoh H. Numerical simulation of the collapse of a bidispersed granular column using DEM and elastoplastic SPH. Comput Part Mech. 2025;12(3):1717–28. doi:10.1007/s40571-024-00896-8. [Google Scholar] [CrossRef]

12. Andersen S, Andersen L. Modelling of landslides with the material-point method. Comput Geosci. 2010;14(1):137–47. doi:10.1007/s10596-009-9137-y. [Google Scholar] [CrossRef]

13. Pastor M, Haddad B, Sorbino G, Cuomo S, Drempetic V. A depth-integrated, coupled SPH model for flow-like landslides and related phenomena. Int J Numer Anal Meth Geomech. 2009;33(2):143–72. doi:10.1002/nag.705. [Google Scholar] [CrossRef]

14. Iverson RM. The physics of debris flows. Rev Geophys. 1997;35(3):245–96. doi:10.1029/97rg00426. [Google Scholar] [CrossRef]

15. Gray JMNT, Wieland M, Hutter K. Gravity-driven free surface flow of granular avalanches over complex basal topography. Philos Trans R Soc A Math Phys Eng Sci. 1999;357(1985):1841–74. doi:10.1098/rspa.1999.0383. [Google Scholar] [CrossRef]

16. Godunov SK. Finite difference method for the computation of discontinuous solutions of the equations of fluids dynamics. Mat Sb. 1959;47(7):271–306. doi:10.7868/s3034503025070108. [Google Scholar] [CrossRef]

17. Harten A, Engquist B, Osher S, Chakravarthy SR. Uniformly high order accurate essentially non-oscillatory schemes III. J Comput Phys. 1987;71(2):231–303. doi:10.1016/0021-9991(87)90031-3. [Google Scholar] [CrossRef]

18. Pudasaini SP, Wang Y, Hutter K. Modelling debris flows down general channels. Nat Hazards Earth Syst Sci. 2005;5(6):799–819. doi:10.5194/nhess-5-799-2005. [Google Scholar] [CrossRef]

19. Lax PD, Richtmyer RD. Survey of the stability of linear finite difference equations. Commun Pure Appl Math. 1956;9(2):267–93. doi:10.1002/cpa.3160090206. [Google Scholar] [CrossRef]

20. Cockburn B, Shu CW. The runge-kutta discontinuous Galerkin method for conservation laws V multidimensional systems. J Comput Phys. 1998;141(2):199–224. doi:10.1006/jcph.1998.5892. [Google Scholar] [CrossRef]

21. LeVeque RJ. Finite volume methods for hyperbolic problems. Cambridge, UK: Cambridge University Press; 2002. doi:10.1017/cbo9780511791253. [Google Scholar] [CrossRef]

22. Bozorgpour R, Darian HM. Recent advancements in fluid flow simulation using the WENO scheme: a comprehensive review. J Nonlinear Math Phys. 2025;32(1):14. doi:10.1007/s44198-025-00269-6. [Google Scholar] [CrossRef]

23. Baeza A, Bürger R, Mulet P, Zorío D. An efficient third-order WENO scheme with unconditionally optimal accuracy. SIAM J Sci Comput. 2020;42(2):A1028–51. doi:10.1137/19m1260396. [Google Scholar] [CrossRef]

24. MiDi G. On dense granular flows. Eur Phys J E. 2004;14(4):341–65. doi:10.1140/epje/i2003-10153-0. [Google Scholar] [PubMed] [CrossRef]

25. Jop P, Forterre Y, Pouliquen O. A constitutive law for dense granular flows. Nature. 2006;441(7094):727–30. doi:10.1038/nature04801. [Google Scholar] [PubMed] [CrossRef]

26. Kamrin K, Hill KM, Goldman DI, Andrade JE. Advances in modeling dense granular media. Annu Rev Fluid Mech. 2024;56(1):215–40. doi:10.1146/annurev-fluid-121021-022045. [Google Scholar] [CrossRef]

27. Johnson PC, Jackson R. Frictional-collisional constitutive relations for granular materials, with application to plane shearing. J Fluid Mech. 1987;176:67–93. doi:10.1017/s0022112087000570. [Google Scholar] [CrossRef]

28. Chassagne R, Bonamy C, Chauchat J. A frictional-collisional model for bedload transport based on kinetic theory of granular flows: discrete and continuum approaches. J Fluid Mech. 2023;964:A27. doi:10.1017/jfm.2023.335. [Google Scholar] [CrossRef]

29. Yuan L, Liu W, Zhai J, Wu SF, Patra AK, Pitman EB. Refinement on non-hydrostatic shallow granular flow model in a global Cartesian coordinate system. Comput Geosci. 2018;22(1):87–106. doi:10.1007/s10596-017-9672-x. [Google Scholar] [CrossRef]

30. Moukalled F, Mangani L, Darwish M. The finite volume method. In: The finite volume method in computational fluid dynamics. Berlin/Heidelberg, Germany: Springer; 2016. p. 103–35. [Google Scholar]

31. Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc. 1977;181(3):375–89. doi:10.1093/mnras/181.3.375. [Google Scholar] [CrossRef]

32. Lucy LB. A numerical approach to the testing of the fission hypothesis. Astron J. 1977;82:1013. doi:10.1086/112164. [Google Scholar] [CrossRef]

33. Yang S, Wang X, Liu Q, Pan M. Numerical simulation of fast granular flow facing obstacles on steep terrains. J Fluids Struct. 2020;99(8):103162. doi:10.1016/j.jfluidstructs.2020.103162. [Google Scholar] [CrossRef]

34. Huang C, Hu C, An Y, Shi C, Feng C, Wang H, et al. Numerical simulation of the large-scale Huangtian (China) landslide-generated impulse waves by a GPU-accelerated three-dimensional soil–water coupled SPH model. Water Resour Res. 2023;59(6):e2022WR034157. doi:10.1029/2022WR034157. [Google Scholar] [CrossRef]

35. Ge W, Li J. Simulation of particle-fluid systems with macro-scale pseudo-particle modeling. Powder Technol. 2003;137(1–2):99–108. doi:10.1016/j.powtec.2003.08.034. [Google Scholar] [CrossRef]

36. Oger G, Doring M, Alessandrini B, Ferrant P. An improved SPH method: towards higher order convergence. J Comput Phys. 2007;225(2):1472–92. doi:10.1016/j.jcp.2007.01.039. [Google Scholar] [CrossRef]

37. Xu F, Zheng MJ, Kikuchi M. Constant consistency kernel function and its formulation. Chin J Comput Mech. 2008;25(1):48–53. (In Chinese). [Google Scholar]

38. Zhang S, Wu D, Lourenço SDN, Hu X. A generalized non-hourglass updated Lagrangian formulation for SPH solid dynamics. Comput Meth Appl Mech Eng. 2025;440(3):117948. doi:10.1016/j.cma.2025.117948. [Google Scholar] [CrossRef]

39. Khan PM, Nayak SK, Nair P. A motif based wall boundary model for smoothed particle hydrodynamics. Phys Fluids. 2025;37(5):053317. doi:10.1063/5.0269343. [Google Scholar] [CrossRef]

40. Monaghan JJ. Simulating free surface flows with SPH. J Comput Phys. 1994;110(2):399–406. doi:10.1006/jcph.1994.1034. [Google Scholar] [CrossRef]

41. Bagheri M, Mohammadi M, Riazi M. A review of smoothed particle hydrodynamics. Comput Part Mech. 2024;11(3):1163–219. doi:10.1007/s40571-023-00679-7. [Google Scholar] [CrossRef]

42. Tahmasebi P. A state-of-the-art review of experimental and computational studies of granular materials: properties, advances, challenges, and future directions. Prog Mater Sci. 2023;138(111):101157. doi:10.1016/j.pmatsci.2023.101157. [Google Scholar] [CrossRef]

43. Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique. 1979;29(1):47–65. doi:10.1680/geot.1979.29.1.47. [Google Scholar] [CrossRef]

44. Zhou ZY, Kuang SB, Chu KW, Yu AB. Discrete particle simulation of particle-fluid flow: model formulations and their applicability. J Fluid Mech. 2010;661:482–510. doi:10.1017/s002211201000306x. [Google Scholar] [CrossRef]

45. Lacaze L, Kerswell RR. Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity. Phys Rev Lett. 2009;102(10):108305. doi:10.1103/PhysRevLett.102.108305. [Google Scholar] [PubMed] [CrossRef]

46. Scaringi G, Fan X, Xu Q, Liu C, Ouyang C, Domènech G, et al. Some considerations on the use of numerical methods to simulate past landslides and possible new failures: the case of the recent Xinmo landslide (Sichuan, China). Landslides. 2018;15(7):1359–75. doi:10.1007/s10346-018-0953-9. [Google Scholar] [CrossRef]

47. Madan N, Rojek J, Nosewicz S. Convergence and stability analysis of the deformable discrete element method. Numerical Meth Eng. 2019;118(6):320–44. doi:10.1002/nme.6014. [Google Scholar] [CrossRef]

48. Kruggel-Emden H, Sturm M, Wirtz S, Scherer V. Selection of an appropriate time integration scheme for the discrete element method (DEM). Comput Chem Eng. 2008;32(10):2263–79. doi:10.1016/j.compchemeng.2007.11.002. [Google Scholar] [CrossRef]

49. Verlet L. Computer “Experiments” on classical fluids. I. Thermodynamical properties of Lennard-Jones molecules. Phys Rev. 1967;159(1):98–103. doi:10.1103/physrev.159.98. [Google Scholar] [CrossRef]

50. Computational Fluid Dynamics Committee. Guide for the verification and validation of computational fluid dynamics simulations (AIAA G-077-1998(2002)). Reston, VA, USA: AIAA Aerospace Research Central; 1998. doi:10.2514/4.472855.001. [Google Scholar] [CrossRef]

51. Berger R, Kloss C, Kohlmeyer A, Pirker S. Hybrid parallelization of the LIGGGHTS open-source DEM code. Powder Technol. 2015;278(2):234–47. doi:10.1016/j.powtec.2015.03.019. [Google Scholar] [CrossRef]

52. Huang C, Sun Y, An Y, Shi C, Feng C, Liu Q, et al. Three-dimensional simulations of large-scale long run-out landslides with a GPU-accelerated elasto-plastic SPH model. Eng Anal Bound Elem. 2022;145(3):132–48. doi:10.1016/j.enganabound.2022.09.018. [Google Scholar] [CrossRef]

53. Yang S, Zhang H, Yu W, Cheng P, Liu Q. Numerical study of interaction between granular flow and an array of obstacles by a bed-fitted depth-averaged model. Chin J Theor Appl Mech. 2021;53:3399–412. doi:10.6052/0459-1879-21-200. [Google Scholar] [CrossRef]

54. Lajeunesse E, Monnier JB, Homsy GM. Granular slumping on a horizontal surface. Phys Fluids. 2005;17(10):103302. doi:10.1063/1.2087687. [Google Scholar] [CrossRef]

55. Gray JMNT, Cui X. Weak, strong and detached oblique shocks in gravity-driven granular free-surface flows. J Fluid Mech. 2007;579:113–36. doi:10.1017/s0022112007004843. [Google Scholar] [CrossRef]

56. Juez C, Caviedes-Voullième D, Murillo J, García-Navarro P. 2D dry granular free-surface transient flow over complex topography with obstacles. Part II: numerical predictions of fluid structures and benchmarking. Comput Geosci. 2014;73(2):142–63. doi:10.1016/j.cageo.2014.09.010. [Google Scholar] [CrossRef]

57. Yin Y, Xing A. Aerodynamic modeling of the Yigong gigantic rock slide-debris avalanche, Tibet, China. Bull Eng Geol Environ. 2012;71(1):149–60. doi:10.1007/s10064-011-0348-9. [Google Scholar] [CrossRef]

58. Wang Z, Zhang C, Haidn OJ, Adams NA, Hu X. FVM realization in Eulerian SPH: a comparative study within unified codebase. Comput Phys Commun. 2025;314(3):109683. doi:10.1016/j.cpc.2025.109683. [Google Scholar] [CrossRef]