A stabilized total-Lagrangian SPH method for large deformation and failure in geomaterials

07/13/2019 ∙ by Md Rushdie Ibne Islam, et al. ∙ BOKU Vienna 0

Conventional smoothed particle hydrodynamics based on Eulerian kernels (CESPH) is widely-used in large deformation analysis in geomaterials. Despite being popular, it suffers from tensile instability and rank-deficiency; thus, it needs several numerical treatments to be stable. In this work, we present a stabilized total-Lagrangian SPH method (TLSPH), which is inherently free of tensile instability. A stiffness-based hourglass control algorithm is employed to cure the hourglass mode caused by rank-deficiency. Periodic update of reference configuration is used in simulations to allow TLSPH to model large deformation and post-failure flow in geomaterials. Several numerical examples are presented to show the performance of the stabilized TLSPH method. The comparison between TLSPH and CESPH are discussed. The influences of hourglass control and configuration update are also discussed and shown in the numerical examples. It is found that the presented stabilized TLSPH is robust and can model large deformation and plastic flows in geomaterials. Particularly, the stabilized TLSPH delivers accurate and smooth stress results.



There are no comments yet.


page 10

page 11

page 12

page 13

page 14

page 15

page 16

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In recent years, meshless Lagrangian particle-based methods, e.g. smoothed particle hydrodynamics (SPH) monaghan2012smoothed , reproducing kernel particle method (RKPM) chen2017meshfree , material point method (MPM) bardenhagen2004generalized , become popular in the modeling of geotechnical problems involving large deformation and post-failure material flow. The reason is that these methods completely or partially avoid the use of mesh; thus, they can conveniently model problems with large deformation, strong material distortion, moving boundaries, and separation of material. Among them, SPH is a well-established truly meshless method having wide application in astrophysics, computational fluid dynamics, and solid mechanics monaghan2012smoothed . It was initially used to simulate seepage failure and liquefaction by Maeda et al. maeda2006development and Naili et al. naili20052d . Then Bui and coworkers bui2008lagrangian ; bui2011slope ; bui2013improved presented a comprehensive SPH framework for the modeling of large deformation and post-failure phenomena in geomaterials. Since then, SPH has been widely used in geotechnical applications such as large deformation bui2014novel , granular flow chambon2011numerical ; peng2016unified , soil-structure interaction wang2014frictional ; zhan2019three , among others.

The previous research and applications of SPH in geotechnical engineering usually employ the conventional SPH, where the search of neighboring particles and computation of kernel functions are based on the deformed current configuration. This type of SPH kernel is named Eulerian kernel in the sense that particles can enter and exit its support domain. Consequently, the commonly-used SPH is termed as conventional Eulerian SPH (CESPH) hereafter. When applied to solid mechanics and geomechanics, the CESPH method has a number of deficiencies which reduce its accuracy and stability ganzenmuller2016hourglass . Among them, the tensile instability swegle1995smoothed and hourglass mode dyka1997stress are the main issues restricting the application of CESPH in solid mechanics. Several numerical treatments were proposed to cure these instabilities. For tensile instability, Monaghan monaghan2000sph and Gray et al. gray2001sph developed the artificial pressure/stress method, where repulsive forces are introduced if there is tensile pressure/stress. The repulsive forces prevent particles from approaching too close to each other hence hinder the formation of particle clumps and artificial cracks. For hourglass mode, Dyka et al. dyka1997stress identified that it is caused by rank-deficiency due to insufficient system integration. Accordingly, they proposed to use additional stress integration points to stabilize CESPH computations. Most geotechnical CESPH simulations employ artificial pressure/stress to alleviate the tensile instability bui2008lagrangian ; bui2011slope ; bui2014novel ; peng2019loquat , whereas the hourglass mode attracts less attention.

Later, the root of tensile instability was identified by Belytschko et al. belytschko2000unified to be the employment of Eulerian kernel function and the updated-Lagrangian formulation. Consequently, a total-Lagrangian SPH (TLSPH) was proposed and proved to be free of tensile instability bonet2002alternative ; vignjevic2006sph . As the Lagrangian kernel has nothing to do with the rank-deficiency, TLSPH still suffers from hourglass mode. Consequently, the path to a robust SPH method is to combine TLSPH with stress points or hourglass control algorithm. However, the stress point method drastically increases the computational cost and a clear rule on how to place and move the stress points is missing, which render it less attracting. Recently, an hourglass control algorithm was proposed for TLSPH based on the hourglass control mechanism used by finite elements with reduced integration ganzenmuller2016hourglass ; ganzenmuller2015hourglass . It is proved that with the hourglass control TLSPH has long term stability, free of tensile instability and hourglass mode, has quadratic convergence rate, and gives rise to accurate stress resultsganzenmuller2015hourglass ; zhan2019stabilized . The last one is a desirable feature which is not yet achieved by other SPH variants in solid mechanics.

Although proved to be a reliable numerical method, TLSPH in its original form is unsuitable for large deformation and failure analysis in geomaterials. The total-Lagrangian formulation cannot handle large deformations and plastic flows, in which the physical connections of materials change under distortion, leading to unphysical deformation gradient and negative Jacobian in the total-Lagrangian formulation. In this work, we employ TLSPH to the modeling of large deformation and post-failure flow in geomaterials. To overcome the shortcomings of the original TLSPH, configuration update is employed with clearly defined criterion. The hourglass control algorithm from Ganzenmuüller ganzenmuller2015hourglass is used to cure the hourglass mode. It is shown in the work that the stabilized TLSPH with configuration update can accurately model large deformation in geomaterials. Additionally, the comparison between the CESPH and TLSPH is carried out. The influences of hourglass control and configuration update are also discussed.

2 Governing Equations and Constitutive Model

2.1 Governing equations

For the modeling of large deformation and failure in geomaterials, the governing equations consisting of mass and momentum conservations are written in the following Lagrangian forms in the current configuration


where , and are the density, velocity, and body force per unit mass. and

are the Cauchy stress tensor and divergence/gradient operator, respectively.

denotes the material derivative. Note that the above governing equations are written in respect with the current configuration , where is the coordinate in the initial configuration, and is a mapping between the initial and current configurations.

The Cauchy stress can be decomposed into two parts: deviatoric stress and isotropic hydrostatic pressure



is a identity matrix. The pressure

is defined as .

2.2 Constitutive model

The elasto-plastic Drucker-Prager (DP) model is employed to model the mechanical responses of geomaterials. The DP model has the following yield surface


where is the second invariant of the deviatoric stress , and are two constitutive parameters, which can be related to cohesion and frictional angle as


In this work, the following plastic potential function with a non-associated plastic flow rule is employed


where is a model parameter related to the dilatancy angle


In computations, a predicted elastic stress is first calculated assuming that there is no plastic effect


where is the stress from the previous time step, is the spin tensor, is the deviatoric strain rate component, . , and are the elastic bulk modulus, the shear modulus, and the time step, respectively. Note that the Jaumann stress rate is used to ensure the objectivity in simulations with large deformation and rotation. The predicted stress is then subjected to plastic correction. If the predicted stress is inside the yield surface, it is accepted as the correct stress; otherwise, a plastic correction should be employed. In the plastic correction, the incremental plastic strain is obtained as


where the plastic multiplier increment is calculated as


Consequently, the corrected stress at the current time step is evaluated as


where , , and represents the predicted variables corresponding to the predicted stress .

3 Conventional Sph for Large Deformation Analysis

3.1 Eulerian kernel-based SPH approximation

In SPH, the problem domain is represented by particles carrying field variables and moving with the material. The approximation of field functions and their derivatives are based on SPH kernel functions. In conventional SPH simulations, the kernel functions are defined in the current configuration with support domains of fixed shape. This conventional kernel function is termed as Eulerian kernel because particles can enter and exit a support domain as the material deforms. Note that although the approximation is based on Eulerian kernels, the conventional SPH is still a Lagrangian method because the particles move with the material. In this section, we briefly introduce the conventional Eulerian kernel-based SPH (CESPH) for large deformation analysis.

In CESPH, a field variable is approximated in the integral form


where denotes the SPH integral approximation, is the support domain of the kernel function , where is the distance in the current configuration and is the smoothing length. In this work, the Wendland C2 kernel function is used wendland1995piecewise


where is a normalization factor, and is the normalized distance.

Similarly, the derivative of the function can be approximated as


where denotes the derivative of the kernel function evaluated at .

The discrete form of the SPH approximation equations can be written as the summation of contributions from all SPH particles in the support domain


where and are the values of field variable at particles and respectively. and denote the mass and density of the -th particle. The kernel function and its derivative are simplified as and hereafter.

3.2 Conventional SPH discretization of governing equations

In CESPH, the governing equations (1) and (2) are discretized in the following forms bui2008lagrangian


where and are the artificial viscosity and the artificial pressure respectively. The artificial viscosity is required for stabilizing computations. The following form of monaghan1983shock is used




where and are two parameters controlling the magnitude of the artificial viscosity, is the relative velocity, is the average speed of sound computed as , with being the elastic modulus of the material. The artificial viscosity should provide stabilized solution without over-damping the system.

For CESPH, the tensile instability is a well-known deficiency which severely limits its application. In slope stability analysis, geomaterials usually have moderate cohesion. It is natural to have some computational areas with tensile stress, where tensile instability develops, which leads to low accuracy or even breakdown of simulations. Therefore, in CESPH, additional numerical stabilization is required. In this work, the widely-used artificial pressure is employed to alleviate the tensile instability in CESPH.

The artificial pressure stabilization adds a positive pressure to prevent particle clumping in the area with tensile stress. It has the following form monaghan2000sph ; gray2001sph .


where is a tuning parameter determined through numerical experiments and , being the initial average particle spacing. is taken as if ; otherwise, it is set to zero. is computed in the same way.

Although artificial pressure can alleviate tensile instability in geotechnical simulations bui2008lagrangian ; peng2015sph ; peng2019loquat , it has two main issues. First, the tuning parameter needs to be determined in each simulation. Simulations with a different numerical resolution, boundary conditions, and material properties may have different values of the tuning parameter. This is undesired because it is inconvenient and introduces arbitrariness, reducing the credibility of the method. The second issue for the artificial stabilization is that in some instances, it cannot prevent tensile instability even with very large ; and further increasing leads to the breakdown of simulations. Consequently, CESPH cannot be used to solve this kind of problems.

4 Stabilized Total-Lagrangian Sph With Configuration Update

Due to the tensile instability, CESPH may have difficulties in slope stability analysis. It is found that the origin of the tensile instability in CESPH is the employment of Eulerian kernels, which are defined based on the current configuration belytschko2000unified . On the other hand, the tensile instability is completely avoided if the kernel functions are computed based on the initial configuration (belytschko2002stability, ; rabczuk2004stable, ). The initial configuration-based kernels are termed as total-Lagrangian (TL) kernels, and the SPH method employing TL kernels are called total-Lagrangian SPH (TLSPH). Because it does not have tensile instability, it is a suitable method for slope analysis with cohesion. However, the original form of TLSPH cannot model problems with large material shearing and distortion due to the use of initial configuration-based TL kernels. It also has an hourglass instability caused by rank-deficiency ganzenmuller2015hourglass . In this section, we first briefly introduce the formulations of TLSPH and propose a configuration update scheme and an hourglass control algorithm to enhance the original TLSPH. These two enhancements enable TLSPH to model problems with large material shearing and distortion without hourglass mode.

4.1 TLSPH formulations

Unlike CESPH in which the particle approximation is based on the current particle position , the field variables and its derivative in TLSPH are evaluated using the reference configuration as


where and are the total Lagrangian kernel function and its derivative evaluated at the reference configuration . is the reference/initial density of the -th particle.

In TLSPH, the conservation Eq. (1) and (2) take the following forms de2013total :


where the values defined in the reference configuration, is denoted by a subscript . is the Jacobian and it is computed based on the deformation gradient matrix as . is the first Piola Kirchhoff stress, which is related to the Cauchy stress as


In the above equation, the deformation gradient matrix is computed as



represents the displacement vector, and

is an identity matrix.

In TLSPH, only the momentum Eq. (28) is solved, as the continuity Eq. (27) is satisfied naturally. The discretized form of Eq. (28) in a total Lagrangian manner is as follows


where is the pullback of the artificial viscosity. As mentioned before, the TLSPH formulation is free from any tensile instability; hence, the use of artificial pressure is not required.

For geomaterials, the same DP model is used in TLSPH simulations. Therefore, the strain rate tensor and spin tensor are needed. To this end, we first evaluate the velocity gradient tensor


where is the rate of deformation gradient obtained as


The discretized forms of the deformation gradient matrix and its rate are


With the velocity gradient , the strain rate tensor and spin tensor can be computed as


The integration of the constitutive model in TLSPH is similar to that in CESPH. Once the Cauchy stress is obtained from the constitutive model, it is converted to the first Piola Kirchhoff stress using Eq. (29), which is then used in the momentum conservation equation (31) to compute the acceleration.

4.2 Update of configuration in TLSPH

The TLSPH in its original form is incapable of modeling plastic flows and post-failure behaviors of geomaterials. This is due to the utilization of the initial/reference configuration in the computation. Whenever there is severe deformation of the material, e.g. particle penetration or material separation, the Jacobian becomes negative, leading to unrealistic prediction and breakdown of simulation. The reason behind this is that when the distortion between the reference configuration and the current configuration is too large, the deformation gradient does not remain a well-conditioned one. The ill-conditioned deformation gradient leads to erroneous values of the field variables like stress and strain in the computation.

To solve this problem, Ganzenmüller et al. ganzenmuller2016hourglass and Leroch et al. leroch2016smooth proposed to update the initial/ reference frame periodically to the current configuration


where the reference particle position and the particle volume are updated. This update does not change the stress tensor, or deformation status, thus preserves stress and deformation history. After the update, all the computations are based on the new reference configuration using the same TLSPH formulations until the next update. It is found that with the configuration update problems with large material distortion and plastic flow can be modeled using TLSPH ganzenmuller2016hourglass ; leroch2016smooth .

An essential consideration of the configuration update is how frequently the reference configuration needs to be updated. In ganzenmuller2016hourglass ; leroch2016smooth , the authors proposed to update the reference configuration when the distance between a pair of interacting particles exceeds one half of the smoothing length. This approach was used in the simulation of a tension test of a hyperelastic patch and the scratching in visco-plastic materials using TLSPH. However, they did not provide much details on how to perform the configuration update and why to choose this specific criterion, which is the key to investigate the performance of the TLSPH method with the configuration update.

As the reference configuration needs to be updated whenever the material distortion is too large, a measurement of distortion is required. In this work, we quantify the distortion in the neighborhood of particle as the maximal normalized displacement in its support domain


This definition can be interpreted as the maximal line strain, with respect to the reference configuration, in the support domain of . At each step, the maximal normalized displacement is computed at each particle, then an overall maximal normalized displacement is obtained as


The overall maximal normalized displacement is employed to quantify the material distortion in the simulation. The reference configuration is updated to the current configuration when the following criterion is met


where is a constant defining the acceptable distortion in the simulation, whose appropriate value and influence on results will be discussed in Section 5. It is found that the criterion with the normalized displacement can more realistically consider the distortion in TLSPH than that used in ganzenmuller2016hourglass ; leroch2016smooth , where the initial distances between particles are not taken into account.

4.3 Stabilization with hourglass control

Being a collocation method having approximation and numerical integration at the same set of particles, SPH suffers from rank deficiency, which usually leads to hourglass mode instability. In CESPH, because particles are often in irregular distribution due to large displacement and distortion, the development of hourglass mode is restrained. Therefore, in CESPH, the problem of hourglass mode is usually negligible, so it does not attract attention. However, in TLSPH, because the computation is based on the reference configuration, where particles generally have regular distribution, the hourglass model develops unboundedly and can eventually lead to poor results or even failure of simulation.

Hourglass mode in TLSPH manifests patterns of particle displacements under which no strain is produced. As a result, no stress can be generated to resist this mode of particle displacement, just like the zero-energy mode in the finite element method with reduced integration. Dyka and Ingel dyka1995approach proposed to add additional stress points for strain and stress computations, which can sufficiently suppress hourglass mode. However, the stress point method gives rise to a much higher computational cost. More significantly, there is a lack of detailed investigation on where to place the stress points and how to move them in simulations with large deformation. Consequently, SPH with stress points appears not to find widespread applications.

Because the TLSPH is entirely free from tensile instability, hourglass mode becomes the dominant source of instability in TLSPH simulations. To this end, a stabilization term for hourglass control is imperative for TLSPH. The hourglass control is particularly crucial for long term simulations, as the hourglass mode develops over time belytschko2000unified . In this work, an hourglass control algorithm proposed by Ganzenmüller (ganzenmuller2015hourglass, ; ganzenmuller2016hourglass, ) is employed. This hourglass control algorithm is used in TLSPH simulation of hyperelastic zhan2019stabilized and elastoplastic leroch2016smooth materials, where accurate results and fast convergence are observed. The basic concept of the hourglass control algorithm is to identify the displacements corresponding to the hourglass mode and then enforce additional forces to minimize these unphysical displacements.

As shown in Eq. (30), the deformation gradient tensor for each particle is computed based on the deformation/displacement of its neighboring particles in the support domain. Ideally, the displacement described by this deformation gradient should be a linear function across the support domain. The hourglass mode, therefore, can be identified as the difference between the actual displacements and the displacements corresponding to the deformation gradient. For instance, for a particle, with reference position and deformation gradient , the ideal linear separation between particle and its neighboring particle can be written as


where . However, as mentioned before the actual separation in the current frame is . The difference between these two linear separations is identified as hourglass mode


Due to symmetry for the -th particle, the error of the pair can be obtained as


To control this hourglass mode, a penalty force along between the -th and -th is employed


where is a penalty parameter taken as 50 ganzenmuller2015hourglass ; ganzenmuller2016hourglass , is the average of elastic models i.e. . For a particle the total hourglass control force can be computed based on averaging over the neighboring particles


The hourglass force is added to the momentum Eq. (31). The total acceleration with hourglass reads


As a summary, Algorithm 1 shows the process of computation in one step using the proposed TLSPH with configuration update and hourglass control. For simplicity, the integration method shown in Table 1 is a simple Euler integration. In actual simulations in this work, a second-order predictor-corrector integration is employed. The details of this integrator can be found in Peng et al. peng2019loquat .

Start: known variables are the reference position , the current position , velocity , and Cauchy stress ;
    1. Compute deformation gradient tensor and its rate for each particle;
    2. Compute the first Piola Kirchhoff stress ;
    3. Obtain total acceleration based on the momentum equation and hourglass control;
    4. Update velocity and the current position ;
    5. Calculate and , update Cauchy stress using the DP model;
    6. Check Eq. (41). If necessary, update the reference configuration and the neighbor list;
    7. Add time .
Algorithm 1 The flowchart of the TLSPH with configuration update and hourglass control

5 Numerical examples

The presented stabilized TLSPH method with configuration update is implemented in the GPU-accelerated open-source SPH code LOQUAT peng2019loquat and is employed to model problems with large deformation and failure in this section. First, the collapse of a column of cohesive soil and its post failure movement are modelled using CESPH and TLSPH. Results from these two methods are compared. The influence of the configuration update in TLSPH and the optimal criterion for updating are investigated. The influence of the hourglass control algorithm is also discussed. Next, the TLSPH is used to evaluate the safety factor of a slope. The results are compared with those from CESPH.

Figure 1: Setup of the soil block for large deformation analysis
Density () Young’s module () Friction angle () Cohesion ()
1850 kg/m 1.5 MPa 25 5 KPa
Table 1: Material properties of cohesive soil
Particle spacing () Smoothing length () Time step () Artificial viscosity parameters ()
0.03 m 0.045 m s (2.5, 2.5)
Table 2: Computational parameters for the soil collapse simulation
Figure 2: Conventional SPH and removal of tensile instability with different artificial pressure parameter

5.1 Large deformation and failure of cohesive soil

This test highlights the accuracy and stability of the present TLSPH framework in modeling large deformation and post-failure flows of cohesive soils. As shown in Fig. 1, a rectangular column of cohesive soil subjected to gravity is considered. The density of the soil particles is kg/m. The material parameters and other computational variables are provided in Table 1 and 2. The DP model is used to simulate the soil. The total simulated time is 6 seconds.

5.1.1 Tensile instability

As mentioned before, conventional SPH suffers from tensile instability, which should be alleviated using the artificial pressure method. However, the constant in the artificial pressure method needs to be calibrated for each simulation with different numerical settings and material properties. For instance, when simulating elastic material, Gray et al. gray2001sph found that is the optimal choice. On the other hand, Bui et al. bui2008lagrangian found that produced best results in SPH simulation of geomaterials.

Figure 3: Final deformed shape free from tensile instability with conventional SPH () and TLSPH

In this work, four CESPH simulations with varying are carried out. The results of the deformed soil body and the distribution of accumulated plastic strain are shown in Figure 2. It is obvious that when no treatment is used (), the results suffer from tensile instability in the form of unphysical cracking and particle clumping. The results are improved with and 0.4, but cracking can still be observed. The tensile instability is completely removed in the simulation with . Further increasing leads to the breakdown of the simulation. It is found that with tensile instability, CESPH could not correctly model the failure pattern, which features a branching in the shear band. This branching of the shear band can only be observed if tensile instability is sufficiently suppressed, as observed in the results from .

Figure 4: Comparison of the stabilised deformed shape with conventional SPH () and TLSPH

The same case is simulated using the proposed TLSPH with configuration update. The comparison between the results from CESPH and TLSPH is shown in Figure 5, where the artificial pressure constant for CESPH is , and the configuration update criterion for TLSPH is . It is observed that with the configuration update, the proposed TLSPH can model problems with large deformation and post-failure flow. It does not need tensile instability treatment and overcomes the difficulty of the original TLSPH in large deformation modeling. Moreover, the presented TLSPH produces results with two distinct shear bands, giving more clear modeling of the progressive failure of the soil column. On the other hand, the shear band in the CESPH results is smeared out on a broader area, indicating that CESPH cannot give detailed modeling of the failure process. This is particularly problematic if one wants to investigate the failure mechanism in more complex cases. The final profiles obtained from the two methods are shown in Figure 4. It is found that the final deformed shapes of the soil body from the two simulations are similar but not without differences. The TLSPH results show a clear step-wise top surface due to the two shear bands. On the other hand, the surface obtained from the CESPH simulation has unphysical bumps caused by artificial pressure.

The deforming and failing process of the soil body obtained using the CESPH and TLSPH are shown in Figure 5. A shear band first initiates at the toe of the column and then propagates into the soil. In the results from CESPH, the shear band is smeared into a wider area, and a branching results in a secondary shear band near the top surface. On the other hand, the second shear band is formed at an early stage in the simulation using TLSPH. Two distinct shear bands are obtained in the TLSPH results. It is found that the TLSPH with configuration update is free of tensile instability and can model large deformation and failure of geomaterials sufficiently. The CESPH, although giving rise to stable simulations with tuned tensile instability coefficients, cannot accurately capture the complete failure process.

Figure 5: Comparison of the deformed shape and accumulated plastic strain at different time step with conventional SPH () and TLSPH

5.1.2 Update of configurations in TLSPH

As mentioned in Section 4.2, the TLSPH in its original form cannot model large plastic deformation, which is common in geomaterials. For instance, Fig. 6 shows the results of the same cohesive soil collapse using TLSPH without configuration update. The simulation terminates at second. This is because the current configuration is much distorted compared with the reference configuration, resulting in unphysical deformation gradient and negative Jacobian at some particles. Therefore, it is necessary to update the reference configuration after a certain degree of distortion. In this work, the overall maximal normalized displacement is employed to measure the distortion. If is larger than a threshold , the reference configuration is updated to the current configuration. That is, the update is performed if the current inter-particle distance is times larger than the original inter-particle distance for at least one pair of particles in the computation.

Figure 6: Effect of update of configurations in TLSPH

To investigate the influence of the threshold , three simulations using TLSPH are performed with , 4, and 8. The results are also given in Fig. 6. With , the simulation terminates at 1.1 s, indicating that the frequency of the configuration update is too low to prevent ill-conditioned deformation gradient. The simulation with is stable, as observed in Fig. 6. However, instabilities can be found in the plastic zones, which means that the deformation gradients are still ill-conditioned. Although this ill-condition does not directly lead to termination of the simulation, it reduces the computational accuracy and results in unphysical deformation patterns. With , the simulation is stable and gives rise to well-defined shear bands. No instability is observed. Therefore, is employed in the following simulations in this work.

It should be pointed out that the choice of the threshold can be linked to the physical properties of geomaterials. During shearing of geomaterials, a physical particle is subjected to the frictional contacts and cohesion forces from surrounding particles; thus, the movements of the surrounding particles contribute to the macroscopic variables such as strain and stress. However, as the shear deformation develops, these links cease to exist during shearing, and new particles get in contact. Therefore, in numerical simulation, new reference configuration needs to be established. This is the physical foundation of the configuration update. However, the precise frequency of configuration update depends on many factors such as grain size, deformation characteristics, boundary conditions, and numerical resolution. Hence, it is difficult to determine the threshold quantitatively using physical properties and numerical conditions. Therefore, in this work, is considered as a numerical constant.

5.1.3 Hourglass control

Previous SPH simulation of large deformation in geomaterials usually did not focus on stress accuracy. In the presented TLSPH, smooth stress results can be obtained if the hourglass control is employed, as demonstrated in Fig. 7. It is observed that without hourglass control, the results show typical hourglass instability, having spurious oscillation of stress over adjacent particle layers. The hourglass instability develops over time, eventually leading to completely distorted particle distribution and poor stress results. More significantly, the distorted particle distribution results in incorrect deformation pattern, as shown in Fig. 8. Due to the hourglass instability, the soil body has unphysical settlements, resulting in an erroneous final profile of the deformed soil block. These unphysical settlements are avoided if the hourglass control is activated. It is found that with hourglass control, the simulation is stable even in long term computations. Moreover, the stress results are smooth; no spurious oscillations can be observed even after a long time of computation. As accurate modeling of stress is crucial in certain applications, e.g. soil-structure interactions, it is a desirable feature for the presented TLSPH method to provide satisfactory stress results.

Figure 7: Comparison of the effect of hourglass control force at time step with TLSPH
Figure 8: Comparison of the stabilised deformed shape in TLSPH with and without hourglass control

5.2 Safety factor analysis

In this section, the TLSPH framework is utilised to evaluate the factor of safety of a homogeneous slope. The strength reduction approach griffiths1999slope is applied to calculate the safety factor. In this approach, the cohesion and friction angle of the soil is reduced incrementally by a factor as


The reduced values of and are then used to check the stability of the slope. The largest value at which the failure of the slope occurs is considered the factor of safety. Once the factor of safety is determined, the post-failure large deformation of the slope is obtained by continuing the simulation.

Figure 9: Setup of the slope for factor of safety analysis
Density () Young’s module () Friction angle () Cohesion ()
1850 kg/m 100.0 MPa 30 5 KPa
Table 3: Material properties of soil for the slope failure set up
Particle spacing () Smoothing length () Time step () Artificial viscosity parameters ()
0.1 m 0.15 m s (2.5, 2.5)
Table 4: Computational parameters for the slope failure simulation

The set up of the slope is shown in Fig. 9, which has the same geometry and boundary condition as in peng2019loquat . The criterion of slope failure is based on the development of the maximum displacement in the slope. This criterion was first proposed by Bui et al. bui2011slope for CESPH analysis of slope stability. Initially, the slope is in an equilibrium state with the original material constants. After the sudden strength reduction, deformations develop in the slope. If the slope remains stable under a reduction factor, the displacement is only caused by the disturbance in material constants. As a result, the maximum displacement should be small and become stable in a short time, and the shape of the time-maximum displacement curve is convex. On the other hand, if the slope cannot remain stable, the maximum displacement increases quickly, and the curve is concave. In this work, the factor is initially taken as and then increased by for each simulation. The maximum displacement for each simulation is recorded over time. If under a certain factor, the time-maximum displacement curve is concave, this factor is taken as the safety factor of the slope.

Figure 10: Distribution of initial stress
(a) Conventional SPH
Figure 11: Displacement over time at different factor of safety

The material constants and the numerical settings are given in Table 3 and 4, respectively. Initially, the original material constants are employed in the computation to get the initial geostatic stresses. The numerical damping approach introduced by Bui et al. bui2013improved is used to accelerate the consolidation process. Once the steady state is achieved, the original particle positions are restored, and all variable except stresses are set as zero. Then the reduced material constants are employed in the simulation to check the stability of the slope. The obtained initial stresses in the slope are given in Fig. 10. It is found using the TLSPH with hourglass control, the accuracy of the stress fields is satisfactory.

Figure 12: Displacement of the slope at different time steps with
Figure 13: Plastic strain of the slope at different time steps with

Figure 11(b) gives the time history of the maximum displacement under varying reduction factor. It is found that if the reduction factor is less than 1.7, the maximum displacement quickly becomes stable and remains a very small value in the whole computation. Correspondingly, the shape of the time-maximum displacement curve is convex, indicating that the slope is stable if the reduction factor is less than 1.7. On the other hand, the maximum displacement develops unboundedly in the interested time span and reaches very large value. The shape of the curve is concave. Consequently, it is found that the slope becomes unstable with reduction factor 1.8. Therefore, for the TLSPH simulation, the obtained safety factor for this slope is 1.8. Additionally, computations using CESPH are performed, and the results are shown in Fig. 11(b). It is observed that the time-maximum displacement curves obtained from the two methods are almost identical. Furthermore, the results of the safety factor are well corroborated by numerical results from literature peng2019loquat .

Figure 14: Stress () distribution of the slope at different time steps with

The developments of displacement and shear band of the slope under reduction factor are given in Fig. 12 and 13, respectively. It is observed that the failure propagates from the toe to the crest of the slope. A well-developed shear band is obtained, whose location and shape are consistent with other numerical results obtained with conventional SPH bui2011slope ; bui2013improved ; peng2019loquat . Above the shear band, a part of the slope body falls and rotates until a new equilibrium state is reached, resulting in concentrated displacements in the sliding body. Although artificial pressure is not used in the TLSPH simulation, no tensile instability can be observed. Furthermore, the distribution of vertical stress at four-time instances are given in Fig. 14. It is found that, except the sliding body, in most part of the slope, the stress does not change. The stress field remains smooth despite the long simulation time. Hourglass instability is avoided, and satisfactory stress results can be obtained, which is a distinct advantage of the presented TLSPH method.

With the configuration update, TLSPH can sufficiently model the post-failure large deformation. Fig. 15 gives the final shape of the slope under reduction factor and 2.0. Reasonably, a higher reduction factor results in larger post-failure deformation. In both cases, large shear deformations are observed along the sliding surface, where particle distortion is high. The original TLSPH is incapable of modeling this problem. Again, in the simulation with , no tensile instability or hourglass mode can be observed, demonstrating the high accuracy and robustness of the presented TLSPH method.

Figure 15: Comparison of the initial configuration (black line) and the deformed slope after s: (a) , (b) .

6 Conclusion

In this work, TLSPH is employed to model large deformation and failure in geomaterials. To improve its accuracy and stability, and to make it be able to simulate problems with large deformation, the following enhancements are used: (1) a stiffness-based hourglass control algorithm is applied to circumvent the hourglass mode caused by rank-deficiency; (2) the reference configuration is periodically updated to the current configuration to prevent numerical instability. A criterion on when to perform configuration update is proposed. An elastoplastic constitutive model is incorporated into the stabilized TLSPH to model geomaterials.

Two cases with cohesive soils involving large deformation and plastic flow are modeled using the stabilized TLSPH method. The results from CESPH and TLSPH show that both methods result in similar deformation patterns provided that appropriate artificial pressure/stress is used in CESPH. However, it is found that TLSPH is free of tensile instability, and it can give more distinct shear bands and more detailed deformation pattern. The criterion on configuration update is discussed, and the optimal choice of the criterion is found to be . The hourglass control is found to be essential for the numerical stability in long-term simulations. With the hourglass control, the stabilized TLSPH can give very accurate and smoothed stress fields, which is otherwise difficult for other SPH variants for solid mechanics. The presented stabilized TLSPH proves to be a sufficient numerical tool for large deformation analysis in geomaterials.


This work receives funding from the European Union’s Horizon 2020 programme under grant agreement No. 778627 and the Austrian Research Promotion Agency (FFG) under the project No. 865963.


  • (1) J. J. Monaghan, Smoothed particle hydrodynamics and its diverse applications, Annual Review of Fluid Mechanics 44 (2012) 323–346.
  • (2) J. S. Chen, M. Hillman, S. W. Chi, Meshfree methods: progress made after 20 years, Journal of Engineering Mechanics 143 (4) (2017) 04017001.
  • (3)

    S. G. Bardenhagen, E. M. Kober, The generalized interpolation material point method, Computer Modeling in Engineering and Sciences 5 (6) (2004) 477–496.

  • (4) K. Maeda, H. Sakai, M. Sakai, Development of seepage failure analysis method of ground with smoothed particle hydrodynamics, JSCE Structural Engineering/Earthquake Engineering 23 (2) (2006) 307s–319s.
  • (5) M. Naili, T. Matsushima, Y. Yamada, A 2D smoothed particle hydrodynamics method for liquefaction induced lateral spreading analysis, JSCE Journal of Applied Mechanics 8 (2005) 591–599.
  • (6) H. H. Bui, R. Fukagawa, K. Sako, S. Ohno, Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geomaterial using elastic–plastic soil constitutive model, International Journal for Numerical and Analytical Methods in Geomechanics 32 (12) (2008) 1537–1570.
  • (7) H. H. Bui, R. Fukagawa, K. Sako, J. C. Wells, Slope stability analysis and discontinuous slope failure simulation by elasto-plastic smoothed particle hydrodynamics (SPH), Géotechnique 61 (7) (2011) 565–574.
  • (8) H. H. Bui, R. Fukagawa, An improved SPH method for saturated soils and its application to investigate the mechanisms of embankment failure: Case of hydrostatic pore-water pressure, International Journal for Numerical and Analytical Methods in Geomechanics 37 (1) (2013) 31–50.
  • (9) H. H. Bui, J. K. Kodikara, A. Bouazza, A. Haque, P. G. Ranjith, A novel computational approach for large deformation and post-failure analyses of segmental retaining wall systems, International Journal for Numerical and Analytical Methods in Geomechanics 38 (13) (2014) 1321–1340.
  • (10) G. Chambon, R. Bouvarel, D. Laigle, M. Naaim, Numerical simulations of granular free-surface flows using smoothed particle hydrodynamics, Journal of Non-Newtonian Fluid Mechanics 166 (12-13) (2011) 698–712.
  • (11) C. Peng, X. G. Guo, W. Wu, Y. Q. Wang, Unified modelling of granular media with smoothed particle hydrodynamics, Acta Geotechnica 11 (6) (2016) 1231–1247.
  • (12) J. Wang, D. Chan, Frictional contact algorithms in SPH for the simulation of soil–structure interaction, International Journal for Numerical and Analytical Methods in Geomechanics 38 (7) (2014) 747–770.
  • (13) L. Zhan, C. Peng, B. Y. Zhang, W. Wu, Three-dimensional modeling of granular flow impact on rigid and deformable structures, Computers and Geotechnics 112 (2019) 257–271.
  • (14) G. C. Ganzenmüller, M. Sauer, M. May, S. Hiermaier, Hourglass control for smooth particle hydrodynamics removes tensile and rank-deficiency instabilities, The European Physical Journal Special Topics 225 (2) (2016) 385–395.
  • (15) J. Swegle, D. Hicks, S. Attaway, Smoothed particle hydrodynamics stability analysis, Journal of Computational Physics 116 (1) (1995) 123–134.
  • (16) C. T. Dyka, P. W. Randles, R. P. Ingel, Stress points for tension instability in SPH, International Journal for Numerical Methods in Engineering 40 (13) (1997) 2325–2341.
  • (17) J. J. Monaghan, SPH without a tensile instability, Journal of Computational Physics 159 (2) (2000) 290–311.
  • (18) J. P. Gray, J. J. Monaghan, R. P. Swift, SPH elastic dynamics, Computer Methods in Applied Mechanics and Engineering 190 (49-50) (2001) 6641–6662.
  • (19) C. Peng, S. Wang, W. Wu, H. S. Yu, C. Wang, J. Y. Chen, LOQUAT: an open-source GPU-accelerated SPH solver for geotechnical modeling, Acta Geotechnica (2019) 1–19.
  • (20) T. Belytschko, Y. Guo, W. Kam Liu, S. Ping Xiao, A unified stability analysis of meshless particle methods, International Journal for Numerical Methods in Engineering 48 (9) (2000) 1359–1400.
  • (21) J. Bonet, S. Kulasegaram, Alternative total Lagrangian formulations for corrected smooth particle hydrodynamics (CSPH) methods in large strain dynamic problems, Revue Européenne des Éléments Finis 11 (7-8) (2002) 893–912.
  • (22) R. Vignjevic, J. R. Reveles, J. Campbell, SPH in a total Lagrangian formalism, Computer Modeling in Engineering & Sciences 14 (3) (2006) 181.
  • (23) G. C. Ganzenmüller, An hourglass control algorithm for Lagrangian smooth particle hydrodynamics, Computer Methods in Applied Mechanics and Engineering 286 (2015) 87–106.
  • (24) L. Zhan, C. Peng, B. Zhang, W. Wu, A stabilized TL–WC SPH approach with GPU acceleration for three–dimensional fluid–structure interaction, Journal of Fluids and Structures 86 (2019) 329–353.
  • (25) H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Advances in computational Mathematics 4 (1) (1995) 389–396.
  • (26) J. Monaghan, R. A. Gingold, Shock simulation by the particle method sph, Journal of Computational Physics 52 (2) (1983) 374–389.
  • (27) C. Peng, W. Wu, H. S. Yu, C. Wang, A SPH approach for large deformation analysis with hypoplastic constitutive model, Acta Geotechnica 10 (6) (2015) 703–717.
  • (28) T. Belytschko, S. Xiao, Stability analysis of particle methods with corrected derivatives, Computers & Mathematics with Applications 43 (3-5) (2002) 329–350.
  • (29) T. Rabczuk, T. Belytschko, S. Xiao, Stable particle methods based on Lagrangian kernels, Computer Methods in Applied Mechanics and Engineering 193 (12-14) (2004) 1035–1063.
  • (30) T. De Vuyst, R. Vignjevic, Total lagrangian sph modelling of necking and fracture in electromagnetically driven rings, International Journal of Fracture 180 (1) (2013) 53–70.
  • (31) S. Leroch, M. Varga, S. Eder, A. Vernes, M. R. Ripoll, G. Ganzenmüller, Smooth particle hydrodynamics simulation of damage induced by a spherical indenter scratching a viscoplastic material, International Journal of Solids and Structures 81 (2016) 188–202.
  • (32) C. Dyka, R. Ingel, An approach for tension instability in smoothed particle hydrodynamics (sph), Computers & structures 57 (4) (1995) 573–580.
  • (33) D. V. Griffiths, P. A. Lane, Slope stability analysis by finite elements, Geotechnique 49 (3) (1999) 387–403.