Abstract:

This paper proposes a novel smoothed particle hydrodynamics (SPH) model to improve the accuracy and stability of large deformation analysis of granular materials. To eliminate the spurious stress field commonly observed in post-failure processes in geotechnical SPH applications, an innovative stress diffusion term is constructed in this study.A novel boundary treatment method is proposed based on a renormalization technique, providing wall boundary conditions with first-order consistency for geotechnical applications; this method can effectively handle complex geometries such as sharp corners.A series of numerical examples verify that: this boundary formulation improves computational accuracy even with relatively coarse particle resolution;the stress diffusion term can significantly reduce the numerical noise affecting the stress field under large deformations, yielding smooth stress distributions.Finally, the developed SPH model is applied to simulations of granular flow and tunnel face collapse, and the results show good agreement with experimental data and demonstrate convergence.

Introduction:

Geophysical flows involving large deformations are closely related to many natural phenomena and engineering applications, including avalanches and rockfalls, debris flows in hazard mitigation, landslides in slope stability, soil liquefaction in seismic design, and internal erosion in dam maintenance. Numerical analysis is a powerful tool for studying progressive failure in geophysical flows and predicting their potentially catastrophic consequences; however, the presence of free surfaces and extremely large deformations are the main challenges in modeling such flows.

In computational geomechanics, the most widely used predictive tools are based on the finite element and finite difference methods. These methods are very efficient in solving partial differential equations (PDEs) and have been widely used in the study of geotechnical engineering problems. Despite the success of mesh-based methods, the use of a mesh still introduces numerical difficulties, such as mesh distortion and domain discretization, which hinder their application to large deformation and free surface problems, such as those encountered in granular flows and post-failure analysis of granular materials. Mesh regeneration techniques offer a possible solution to the mesh distortion problem; however, they are computationally expensive and not always practical.

Compared with mesh-based methods, meshless methods provide a more effective solution for large deformation problems, in which mesh connectivity is eliminated. Examples of meshless methods used for granular materials and geotechnical applications include smoothed particle hydrodynamics (SPH), the material point method (MPM) and the particle finite element method (PFEM). Among these, SPH can be regarded as a truly meshless, Lagrangian method.In SPH, the physical domain is discretized into a set of Lagrangian particles that carry macroscopic field properties and move with the material velocity. The interactions between particles are governed by a weight function called the smoothing kernel function.Since local interpolation is used to evaluate the field variables of each particle, no background mesh is required to integrate the governing PDEs. These features make SPH highly suitable for problems involving strong nonlinearity, free surfaces, and large deformations.

Although SPH was originally developed by Gingold and Monaghan and Lucy for astrophysical problems, over the past few decades it has been applied to a wide range of fields, including fluid dynamics, solid fracture, multiphase flow, soil-water interaction, and debris flow. The application of SPH in computational geomechanics began with Bui et al. In their work, an elastic-perfectly plastic Drucker-Prager model with associated and non-associated plastic flow rules was successfully incorporated into the SPH framework to study the post-failure flow of cohesive and non-cohesive soils. Subsequently, various rigorous constitutive models, such as the Bingham flow model, rate-dependent viscoplastic models, hypoplastic models, critical-state-based constitutive models, have been incorporated into the SPH framework to simulate different flow behaviors. In addition, several studies have been devoted to eliminating the shortcomings of SPH in large deformation geomechanical problems, such as tensile instability, by combining SPH with other methods, such as Taylor-Galerkin SPH and Stress-Particle SPH.

In most existing geotechnical SPH applications, the focus has mainly been on kinematic predictions, and satisfactory results have been obtained compared with experimental observations. However, when analyzing stress distributions, high-frequency short-wavelength oscillations are observed under large deformations. Similar to the pressure fluctuations observed in weakly compressible SPH for fluids, this spurious stress field affects the accuracy of the SPH method and reduces its predictive capability, because accurate stress prediction is important for some practical applications (e.g., simulating the impact of debris flows on structures). The stress regularization technique recently proposed by Nguyen et al. provides a solution to the spurious stress field by performing a filter over the entire domain to re-evaluate and redistribute the stress of each particle at a specific temporal frequency.

Nguyen et al. (2017) suggested using moving least squares (MLS) interpolation for stress regularization to achieve first-order correction. In their work, the regularization technique successfully filtered high-frequency numerical noise in the stress field. However, a common problem with using filtering techniques in SPH is that, for long-time simulations, field variables may be over-filtered, thereby affecting the physical mechanisms involved in the simulated process. In addition, the filtering process itself requires an additional particle sweep over the entire domain, which is computationally expensive, especially for large-scale problems. Therefore, it is necessary to seek an alternative approach for stress smoothing that can simultaneously preserve momentum conservation and consistency.

Another important issue associated with SPH is the treatment of boundary conditions.The accuracy and robustness of numerical simulations depend heavily on the performance of the boundary conditions adopted in the SPH model. For SPH applications in geomechanics, the most widely used boundary treatment is the method proposed by Bui et al. (2008). In their work, for soil particles near the boundary, the stress values of neighboring boundary particles are set equal to those of the soil particle, thereby giving a locally uniform stress distribution near the solid boundary. The no-slip boundary condition is imposed by constructing artificial velocities for boundary particles as described in Morris et al. (1997). Subsequent studies have proposed several modifications to the Bui boundary condition. Peng et al. (2015) argued that Bui's treatment of boundary particle stresses may not guarantee the non-penetration condition, and proposed taking the diagonal components of the boundary particle stress tensor as the maximum values of the soil particle of interest, while setting the off-diagonal components equal; this ensures that the boundary generates sufficiently large pressure. Chalk et al. proposed a different modification to the same problem, employing a single layer of dummy wall particles to apply repulsive forces to ensure the non-penetration condition. Another possible boundary treatment relies on the use of dummy particles and renormalization techniques. The field quantities (e.g., stress, velocity) of dummy boundary particles are extrapolated from their neighboring soil particles through a Shepard-corrected kernel function to give zeroth-order consistency. It is also worth mentioning the work of Douillet-Grellier et al. and Zhao et al. on developing stress boundary conditions without boundary particles. The basic idea of their treatment is to complete the truncated kernel support of particles near the boundary by including an additional stress term in the momentum equation. This type of Dirichlet wall boundary condition can be applied when boundary stresses are required, for example, in simulations of triaxial tests, Brazilian tests, and coin-shaped crack problems.

The boundary treatment techniques available for SPH geotechnical applications can ensure at most zeroth-order consistency. The lack of consistency in boundary treatment may lead to poor predictions of the stress field near boundaries. In addition, although many of these boundary methods perform well in test cases with simple geometries, their performance has not been evaluated under more complex geometries.

This work develops a new approach that can minimize the spurious stress field commonly encountered in SPH simulations of geomechanical problems and address some limitations of existing boundary conditions. Inspired by the density diffusion scheme in fluid dynamics, which provides a consistent numerical algorithm to eliminate pressure noise and improve stability, this study considers a diffusion process to improve the stress field.To accurately resolve the stress distribution near boundaries while being able to handle complex geometries, this study modifies the fixed ghost particle technique proposed by Marrone et al. to treat wall boundary conditions. To the best of the authors' knowledge, this is the first time such a boundary technique has been extended to simulate elastoplastic solids with stress in SPH. Unlike Marrone et al., who used MLS interpolation to evaluate flow quantities, this paper adopts the corrected SPH interpolation proposed by Liu and Liu (2006) to restore first-order consistency at the boundary. This process can be achieved within a single loop, thereby minimizing computational cost.

The remainder of this paper is organized as follows: Section 2 briefly describes the governing equations and constitutive model. Then, Section 3 presents the SPH discretization of the governing equations and the development of the numerical diffusion term. Section 4 describes in detail the first-order consistent wall boundary treatment. In Section 5, the accuracy and robustness of the new SPH scheme are evaluated through several verification cases, including Couette flow, simple shear test, static soil with a wedge, two-dimensional granular column collapse, and tunnel face collapse. Finally, Section 6 presents the main findings and conclusions.

Conclusions:

This paper proposes a novel SPH method for large deformation analysis of geomaterials. A novel no-slip/free-slip boundary treatment method is introduced that employs the fixed ghost particle technique combined with first-order consistent interpolation.The application of such a boundary method in SPH for solving geotechnical engineering problems is the first of its kind. This new boundary formulation is suitable for complex geometries and is shown to provide accurate stress and velocity predictions. To address the spurious stress field under large deformations, a common problem in SPH simulation of geomaterials, an improved scheme is proposed by introducing a numerical diffusion term into the stress. Several methods for constructing the diffusion term are discussed in this paper and compared with the new diffusion term algorithm.The results show that the algorithm can smooth numerical noise while maintaining consistency at the boundary. This treatment is easy to implement and does not require an additional, costly stress interpolation process.Simulation results of Couette flow and static dry soil in complex geometries show convergent and accurate predictions of velocity and stress distributions, demonstrating the accuracy, robustness, and applicability of the proposed technique. The new scheme is finally applied to simulations of granular column collapse and tunnel face collapse, and the results show good agreement with experimental results and provide smooth stress fields.