Abstract:

A new method combining Smoothed Particle Hydrodynamics (SPH) with Finite Element Analysis (FEA) is proposed to address the computational challenges of simulating large deformation processes in the interaction between granular media and rigid solids.The introduction of ghost particles at free surfaces, solid boundaries, and the SPH/finite element interface preserves particle sampling accuracy near boundaries and interfaces.The proposed method of applying an initial stress field in the SPH/FEA model enables simulation of loading tests under confined conditions, which has rarely been simulated with SPH in previous studies.The SPH/FEA method is evaluated by simulating biaxial compression tests, shallow foundation bearing capacity, and wheeled locomotion in granular media. The biaxial compression test results are generally consistent with reference solutions established using the finite element method alone. The SPH/FEA simulations yield satisfactory estimates of shallow foundation bearing capacity and appropriately capture the soil failure mechanism of a slope subjected to distributed loading. Because theSPH/FEA method can successfully capture contact or loss of contact between a granular matrix and an intruder (such as a wheel) moving at a prescribed velocity, the dynamic simulation of wheeled locomotion in granular media is more realistic than the velocity field predictions obtained through Resistive Force Theory (RFT)

Introduction:

Manygeotechnical engineeringproblems involve large deformations, and developing adequate computational methods to solve these problems is crucial for engineering practice. In the Discrete Element Method (DEM), the geological material of interest is represented by particle elements. Based on the contact forces generated to ensure equilibrium, the displacement and velocity fields of the particles are computed. Each particle is subjected to gravitational acceleration, elastic contact forces, and dissipative normal and frictional forces from neighboring interacting particles. DEM is well suited for analyzing small-scale soil mechanics problems such as granular column collapse and soil/structure interaction, but its application to high-resolution problems is limited by its computational cost. DEM models require numerous parameters, such as particle size, friction coefficient, and particle shape. Although DEM models are sensitive to these parameters, there is currently no systematic method to efficiently calibrate DEM models. For example, reproducing elastoplastic, hypoplastic, or viscoelastic behavior requires substantial computational effort. Soil/structure interaction models use DEM combined with the finite element method, but the use of DEM/FEM approaches remains hindered by calibration challenges and computational cost.

The Material Point Method (MPM) also employs particle elements to discretize the domain, but unlike DEM, a material point is represented by a set of adjacent particles that are assigned a continuum mechanics constitutive model. The coupled differential equations at the material points can then be solved for field variables such as displacement and stress. MPM has been used to solve complex problems such as saturated clay intrusion or sand column collapse. Although MPM has been successful in studying soil-structure interaction problems, it has certain limitations. First, field variables are defined on a background grid, which functions similarly to the mesh in Finite Element Analysis (FEA). In most MPM software packages, the background grid remains fixed once the simulation begins, which is a limitation for problems involving complex geometries and large deformations. Second, MPM requires modeling the entire soil domain using uniformly sized particles. Third, MPM is known to produce inaccurate results for small-deformation problems in solid mechanics. Using FEM in regions where small deformations are expected can improve computational accuracy in some cases compared with using a pure MPM model. However, coupling the material point model mesh with the finite element model mesh affects the positions of all particles in the grid cells that are in contact with the finite element mesh, which can sometimes produce unwanted particle motion.

Because Smoothed Particle Hydrodynamics (SPH) has a truly meshless nature, such inaccuracies do not arise when FEM is coupled with SPH. Like MPM, SPH also employs particle elements, but the particles are not bound to a grid, which allows dynamic particle refinement and coarsening. SPH is a continuum meshless numerical method that is well suited to large deformation problems involving flowing materials such as fluids or granular media. SPH is based on a Lagrangian description of fluid motion, in which the fluid is represented by a set of particles that move according to the momentum balance equations. Solving these equations requires searching for the nearest neighbors of each particle. DEM is highly accurate in simulating individual particle interactions and resolving detailed contact mechanics and is well suited to small- to medium-scale systems; however, its computational cost scales poorly with particle number, limiting its applicability to large-scale problems. Studies have found that DEM and SPH achieve similar accuracy for classical granular flow problems such as granular column collapse and particle sedimentation. MPM combines Lagrangian and Eulerian approaches and is advantageous for capturing large deformations and complex boundary interactions, but it may suffer from numerical issues such as particle clustering or inaccurate stress transfer. Some features of the MPM and SPH methods are discussed in further detail in Arduino (2018) and Fern et al. (2019). As a fully Lagrangian particle-based method, the SPH method performs best in handling free-surface flows and highly deformable systems, providing a robust solution for continuous granular flow; however, its accuracy may be limited by kernel smoothing artifacts, and its computational cost, although lower than DEM for large-scale systems, is still significant. Efficiency depends not only on the computational method but also on its implementation. Therefore, the choice of software (ANSYS, a MATLAB-based laboratory-developed code, or others) and hardware (number of cores, CPU vs. GPU, etc.) is also very important. In this study, the code leverages the CPU parallelization provided by MATLAB.

Various soil constitutive models have been adopted in SPH formulations to study large-deformation geomechanics problems. SPH has demonstrated remarkable accuracy for dynamic problems in compressible and nearly incompressible media. In particular, Riemann-solver-based spatial discretization can avoid non-physical hydrostatic pressure oscillations, hourglassing, and numerical issues associated with lack of conservation, consistency, stability, or convergence. In addition, a numerical stress diffusion equation effectively eliminates spurious stress oscillations in SPH simulations. Because SPH was originally developed for fluid dynamics, several research groups have used it to simulate soil liquefaction. For example, the pore fluid in saturated water-bearing soil is modeled with SPH, where each fluid particle is subjected to a weighted drag force computed from the displacement field of the solid skeleton, which is itself modeled using the Discrete Element Method (DEM). In another SPH approach, elastoplastic SPH soil particles are used, and particles exceeding a certain plastic strain level are converted into particles obeying non-Newtonian rheology. Biot's poromechanics equations are implemented using a semi-implicit SPH algorithm and allow simulation of slope failure. In this approach, Darcy's drag is updated implicitly, mitigating the permeability-related time-step restriction. SPH has been extended to unsaturated soils by assigning constitutive behavior for effective stress and suction to SPH particles. An SPH model has also been developed to simulate the response of saturated soil to seismic waves, in which flexible foundations and free-field boundaries are specifically used for seismic analysis. Despite recent advances based on parallelization and Graphics Processing Unit (GPU) usage, SPH remains relatively computationally expensive. Moreover, under quasi-static conditions, SPH performs less well than FEM; therefore, SPH is less attractive than FEM for problems with low loading rates relative to the material deformation rate. Indeed, the accuracy of SPH for solid mechanics problems remains an ongoing topic of research and debate, particularly under quasi-static conditions. To apply SPH in solid mechanics, the momentum balance equations are modified to account for solid behavior, which is modeled by constitutive equations for elasticity, plasticity, or damage. By coupling SPH/DDA (Discontinuous Deformation Analysis), it is possible to simulate the interaction between soil (modeled with SPH) and rigid rock blocks (modeled with DDA) during landslides, or to simulate the interaction between fluid and soil volumes undergoing rigid-body motion. Simulating deformable solid surfaces or interfaces between deformable solids and granular media with SPH remains challenging, which limits the soil-structure interaction problems that SPH can investigate.

Herein, we propose a new method that combines the SPH method with FEM and implement an SPH/FEA model in MATLAB to simulate the interaction between a rigid deformable solid (modeled with FEM) and a host granular medium (modeled with SPH) undergoing large deformations under quasi-static conditions.The SPH governing equations are described in Section 2. The principles of the proposed SPH/FEA coupling resolution algorithm are presented in Section 3. Section 4 details the implementation of SPH boundary conditions and initial stress field conditions. We use the proposed SPH/FEA to simulate various stress paths and boundary value problems of interest in geomechanics, such as biaxial compression problems, shallow foundation bearing capacity problems, and wheeled locomotion problems on soil. The results are discussed in Section 5.

Conclusions:

A robust MATLAB solver is proposed that effectively combines Smoothed Particle Hydrodynamics (SPH) and Finite Element Analysis (FEA) for simulating various types of geomechanics problems.The interaction between the SPH and FEA domains is handled by a penalty method, i.e., the particle penetration distance is used to compute the reaction exerted when the FE solid contacts the particles. The penalty method ensures mass conservation.

The introduction of ghost particles at free surfaces, solid boundaries, and the SPH/FEA interface preserves particle sampling accuracy near boundaries and interfaces. Another important contribution of this paper is the proposed method of applying an initial stress field in the SPH/FEA model, which enables simulation of loading tests under confined conditions, which has rarely been simulated with SPH in previous studies.

The effectiveness of the SPH/FEA method is verified by benchmarking against a shallow foundation model based on Terzaghi's bearing capacity equation. The method is then used to simulate intrusion and propulsion. The SPH/FEA method successfully captures the decrease in foundation bearing capacity as the slope is approached, and the predicted bearing capacity factors are consistent with the literature. SPH/FEA also allows simulation of wheeled locomotion on a granular substrate and captures the transition between slip and propulsion as a function of the applied angular velocity at a constant applied horizontal velocity. We consider the SPH/FEA predictions of the soil velocity field to be more realistic than those obtained through Resistive Force Theory (RFT), because SPH/FEA can successfully capture contact or loss of contact between the solid wheel and the granular substrate at different rotational speeds.

The proposed SPH/FEA method leverages the MATLAB Parallel Computing Toolbox. The simulations are run on a 40-core CPU.This setup has acceptable performance for medium-scale problems. The simulation time under this setup scales approximately linearly with the number of SPH particles. This is mainly because the simulation is quasi-static, which avoids performing a full neighbor search at every time step. This strategy significantly reduces the computational cost of the full neighbor search in the SPH domain and the intrusion detection at the SPH/FEA interface. Although this setup is effective for 2D or 3D problems that can be simplified to 2D problems, full 3D simulations will require handling many more particles. The model can potentially be improved by migrating the computation to a GPU-based architecture to exploit its parallel processing capability. Although MATLAB supports GPU computing through gpuArray and a set of GPU-compatible functions, the code developed in this study contains many custom functions that are currently incompatible with GPU execution. In addition, CPU-GPU hybrid approaches in MATLAB often suffer from inefficiencies due to data transfer bottlenecks between the CPU and GPU.Therefore, the SPH/FEA method can be extended to handle the following: (i) migration to full GPU computing in the future as MATLAB gradually enhances its GPU support for complex numerical operations; or (ii) scaling out to high-performance computing clusters with hundreds of CPU cores, which can be readily deployed to large and complex problems.

Although the model proposed in this paper is not benchmarked against other SPH models (such as the SPH/DEM method), doing so in future work could be interesting. The development and validation of the SPH/FEA solver pave the way for its application to more complex engineering problems, such as studies of tool/soil interaction during excavation, modeling of anchor installation and use, and simulation of underground motion of articulated solids. It is worth noting that the MATLAB program developed in this study leverages the advantages of parallel processing and GPU acceleration and is openly available.