Numerical analysis of slope collapse using SPH and the SIMSAND critical state model

Zhao Lu ,Zhuang Jin ,Panagiotis Kotronis

a Civil and Environmental Engineering,University of Macau,Macau,China

b Research and Development Center,Shenzhen Foundation Engineering Co.Ltd.,Shenzhen,518000,China

c Department of Civil and Environmental Engineering,The Hong Kong Polytechnic University,Hong Kong,China

d Southern Marine Science and Engineering Guangdong Laboratory(Guangzhou),Guangzhou,510000,China

e Department of Ocean Science and Engineering,Southern University of Science and Technology,Shenzhen,518000,China

f Institut de Recherche en Génie Civil et Mécanique(GeM),École Centrale de Nantes,Universitéde Nantes,Rue de la Noë,Nantes,44321,France

Keywords:Granular material Smoothed particle hydrodynamics(SPH)Large deformations Landslide Critical state Slope failure Sand

ABSTRACT Geological disasters such as slope failure and landslides can cause loss of life and property.Therefore,reproducing their evolution process is of great importance for risk assessment and mitigation.The recently developed SIMSAND critical state sand model combined with the smoothed particle hydrodynamics(SPH)method is adopted in this work to study slope failure under large deformations.To illustrate the ef ficiency and accuracy of the SIMSAND-SPH approach,a series of slope collapse studies using the discrete element method(DEM)considering various particle shapes(i.e.spherical,tetrahedral and elongated)is adopted as benchmarks.The parameters of the SIMSAND model are calibrated using DEM triaxial tests.In comparison to the DEM simulations,the runout distance and final slope height are well characterized with the SIMSAND-SPH approach with less computational cost.All comparisons show that the SIMSAND-SPH approach is highly ef ficient and accurate,which can be an alternative numerical tool to simulate real scale granular flow.

1.Introduction

Geological disasters accompanied by large deformations and failure of geomaterials and geo-structures occur regularly around the world,e.g.landslides(Liu et al.,1995;Wang et al.,2021),snow avalanche/rockfalls(Bovet et al.,2010),debris flow(Jeong et al.,2015),dam breaks(Shakibaeinia and Jin,2011),soil liquefaction(Di and Sato,2003),seepage damage(Maeda et al.,2006;Yang et al.,2019a),tunnel excavation(Shi et al.,2019)and dynamic erosion(Yang et al.,2019b,2020;Lei et al.,2020).Large-scale solidfailure sliding plays an important role in these geological hazards that generally cause huge destructions and present a long-distance impact.Reproducing their evolution process is therefore necessary for risk assessment and mitigation.

In recent years,the discrete element method(DEM)has been often used to numerically reproduce granular collapse(Lacaze et al.,2008;Soundararajan,2015;Utili et al.,2015)and landslides(Barla et al.,2012;Scholtès and Donzé,2012;Lu et al.,2014).Unlike the finite element method(FEM),DEMtreats the material as an assemblage of rigid discrete particles.Contact forces are used to reproduce the particle interactions and thus the physical nature of the granular body is taken into account.It is worth noting that the number of particles in the DEM simulations should be very small.Therefore,parallel computational strategies are often necessary in order to reduce the computational cost.However,the number of particles in the majority of DEMsimulations is limited for computational ef ficiency and often far from a realistic physical model.In view of this,coupling the DEM with other numerical methods(such as finite elements,finite differences and the material point method(MPM))is regarded as an alternative approach.Taking the FEM-DEMcoupling for example,the granular materials presenting large deformations are modeled by DEM particles,whereas finite elements are used for the domain presenting small deformations.This allows the number of DEMparticles to be greatly reduced,improving the computational ef ficiency,so that the model can be applied to real case studies.Nevertheless,the computational ef ficiency of the DEM is relatively low compared to phenomenological constitutive models,such as the classical Mohr-Coulomb model(Coulomb,1973)or the more advanced Modi fied Cam Clay model(Roscoe and Burland,1968),the SIMSAND(simple sand)model(Yin et al.,2016;Jin et al.,2017a)and the 3Delastoplastic model(Hong et al.,2017,2020).

Another important limitation stems from the fact that severe deformations of the mesh limit the application of the FEM.In contrast,mesh-free methods such as smooth particle hydrodynamics(SPH)and MPM are good alternatives to simulate landslides(McDougall and Hungr,2004;Haddad et al.,2010;Cascini et al.,2014;Llano-Serna et al.,2016;Soga et al.,2016;Peng et al.,2019;Ren et al.,2019;Xiong et al.,2021).For example,SPH has been successfully adopted to analyze the impact of landslideinduced tsunami on shipping in reservoir(Wang et al.,2021).However,numerical simulations considering mesh-free methods and state-dependent advanced constitutive models are still seldom in the literature(Jin et al.,2018a;Yin et al.,2018a).

In order to ef ficiently reproduce slope collapse taking into account some advanced features such as the in fluence of soil density and large deformations,the SIMSAND critical state constitutive model and the SPH method are hereafter coupled.In order to validate the approach,DEM slope collapse simulations are used as benchmarks.For the DEM simulations,three different particle shapes are considered:spherical shape(SS),tetrahedral shape(TS)and elongated shape(ES).The SIMSAND-SPH parameters are calibrated from DEMtriaxial tests using an optimizationbased identi fication method.Comparisons of the DEM and the SIMSAND-SPH simulations illustrate the performance of the approach and its ability to take into account the effect of the particle shape and the soil density on the final deposit con figuration.

2.Numerical model

2.1.Smoothed particle hydrodynamics(SPH)

The SPH method,implemented in the finite element code ABAQUS/Explicit,is adopted due to its ability to simulate large deformation problems.This method was introduced in 1977(Gingold and Monaghan,1977),initially for problems in astrophysics.Further developments extended SPH to solid mechanics where its main advantage is the possibility of dealing with larger local distortions over mesh-based methods.

The core of SPH is interpolation,and all macroscopic variables(such as density,pressure,speed,etc.)are conveniently expressed as integral interpolation calculations by means of a set of values at discrete points.The field variablef(x)at pointxin the numerical domainΩis obtained by the following equation accounting for the effect of neighbouring particles:

whereWis a weighting function,the so-called kernel or smoothing function.Particles in the numerical domain interact with each other through the kernel function with a characteristic radius known as the“smoothing length”,represented byhin Eq.(1).In other words,the physical quantity of any particle is obtained by summing the relevant properties of all the particles that lie within the range of the kernel.The functionf(x)in Eq.(1)is further approximated by the summation over neighbouring particles as

whereVi,ρiandmiare the volume,density and mass of the particlei,respectively;andNis the number of in fluencing particles.The derivatives of the above interpolated integrals are obtained by ordinary derivation without the use of finite differences and without grids,which is the main advantage of SPH over conventional numerical methods(finite differences or finite elements).

By assigning at each particle its own smoothing length and allowing it to vary with time,a SPH simulation is automatically adapted on local conditions.In a very dense region where particles are close,the smoothing length can be relatively short,yielding a high spatial resolution.Conversely,in low density regions where individual particles are far apart and the resolution is low,the smoothing length can be increased,optimizing the computational time.Since the particles can interact and then separate during calculation,SPH is able to deal with very large deformations(Baranowski et al.,2014).

2.2.Explicit finite element method(FEM)

The SPH method implemented in ABAQUS uses the explicit time integration method(Hibbitt et al.,2014).The equilibrium condition is first written with the balance of internal and external forces:

where M is the mass matrix,¨u is the acceleration,P is the applied external force vector,I is the internal force vector,and the subscripttrefers to the time in an explicit dynamic step.

The equations of motion for the body are then integrated using the explicit time central-difference integration rule as follows:

where u is the displacement,˙u is the velocity,andΔtis the time increment.For stability reasons,the time incrementΔtshould be smaller than a limited value,i.e.Δt≈Lmin／cd,whereLminis the smallest finite element dimension andcdis the dilatory wave speed.The incremental displacement is used to calculate the incremental strain.The stresses are updated using the constitutive model and then the internal forces up to a new equilibrium state.

2.3.SIMSAND,a critical state-based model for sand

The SIMSAND model is based on the classical Mohr-Coulomb model with some additional advanced features:the critical state concept(Yin et al.,2016;Jin et al.,2017a),nonlinear elasticity,nonlinear plastic hardening,interlocking effects and a simpli fied threedimensional(3D)strength criterion.The model is able to capture both the state-dependent peak strength and stress-dilatancy(contraction or dilation)(Jin et al.,2017a).The basic constitutive equations are summarized in Table 1.The model parameters and their de finitions are summarized in Table 2.Their calibration can be carried out using either trial-error(Wu et al.,2017)or optimization methods(Jin et al.,2016a,b,2019a).

Table 1Basic constitutive equations of the SIMSAND model.

Table 2Parameters of the SIMSAND model used in slope collapse simulations.

The SIMSAND model was implemented into ABAQUS/Explicit as a user-de fined constitutive model,via the user material subroutine VUMAT.The implementation of the model follows the procedureshown in Fig.1(Hibbitt et al.,2014).The cutting plane algorithm(Ortiz and Simo,1986)is adopted for the stress integration.The implementation was veri fied by simulating drained triaxial tests.

Fig.1.Flowchart of explicit analysis of the SIMSAND model combined with the SPH method in ABAQUS/Explicit.

Fig.2.Schematic plot and spatial discretisation of the slope.

3.Slope collapse simulation using discrete element method(DEM)

In order to validate the combined SIMSAND-SPH approach,a series of DEMsimulations is first conducted and used as benchmark data.More speci fically,the 3D discrete element software PFC3Dcoupled with the finite difference code FLAC3Dare used for the simulations.Fig.2 shows the schematic diagram and spatial discretisation of the slope that has the following dimensions:height is 1 m,top width is 2 m,bottom width is 6 m,and an angle is 45°.The slope thickness in theydirection is 0.1 m,which is about 10 times the average particle size(d50).In the discrete particle region,the total number of particles is about 100,000.The slope is divided in two regions,as shown in Fig.2:the PFC model region where discrete particles are used for the discretisation,and the FLAC model region where finite elements are used for the discretisation.The corresponding areas(PFC and FLAC areas)of the model are built to match the characteristics of slopes in the reality.

Fig.3.Particle shapes:(a)Spherical,(b)tetrahedral,and(c)elongated.

In order to study the effect of particle shape on landslides,three kinds of particle shapes are studied:SS,TS and ES,as displayed in Fig.3.Furthermore,a“dense”slope and a“loose”slope are considered.For the“dense”slope,the initial friction coef ficient of the particles is considered equal to 0.Under the action of gravity,the particles can thus fall and deposit more easily,thereby generating a dense slope.For the“loose”slope,the initial friction coefficient of the particles is taken equal to 1.Table 3 lists the parameters of the 8 studied slope types.In the discrete particle region (PFC model),the contacts between the particles are described with a linear contact model,as shown in Table 4.In the finite volume area(FLAC model),the classical Mohr-Coulomb constitutive model is adopted,calibrated from drainage triaxial tests(Table 5).A gravity field is applied to the entire model.

Table 3Parameters of the discrete particles for the 8 slopes.

Table 4Inter-discrete particle properties.

Table 5Parameters in the finite volume area.

Table 6Calibrated parameters of the SIMSAND constitutive model.

4.Slope collapse simulation using the combined SIMSANDSPH approach

4.1.Numerical model

The geometric size of the DEM calculations is followed in the SIMSAND-SPH approach.The discrete particle and finite volume regions are now replaced by SPH particles and conventional Lagrangian finite elements,respectively.For the SPH domain,a cell size of 0.015 m is estimated by checking mesh-dependency and the total number of particles is 55,300.In the Lagrangian domain,the element size and number are consistent with those of the FLAC3Dmodel.The inner restraining wall is supposed to be smooth and the friction between particles and inner surface of restraining wall is negligible.The bottom of the slope is considered fixed while symmetrical conditions are assumed for the four lateral boundaries.In the finite element area,the classical Mohr-Coulomb constitutive model is again adopted and calibrated from drainage triaxial tests(see Table 5),while the SIMSAND model is considered in the SPH domain.A gravity field is applied to the entire model.

It is worth noting that compared to the DEMcases which require an average simulation time of 24 h for each simulation by adopting 8 CPU cores(processor clock speed equals 2.5 GHz),the average computation time for each SPH case approximately equals 2 h.Moreover,several cases can be simultaneously calculated on the TAIYI server in the Center for Computational Science and Engineering at Southern University of Science and Technology,Shenzhen,China.Therefore,the calculation ef ficiency of proposed method is acceptable.

4.2.Calibration of the SIMSAND model parameters

Fig.4.Comparison between DEM triaxial tests(red symbols)and the SIMSAND constitutive model for different particle shapes(blue lines):(a)SS(rolling resistance=0);(b)TS(rolling resistance=0);(c)ES(rolling resistance=0);and(d)SS(rolling resistance=0.1).

A series of DEM drained triaxial tests is conducted to calibrate the parameters of the SIMSAND constitutive model.For the three kinds of particle shapes,a total of 6 tests,i.e.3 tests on dense samples and 3 tests on loose samples,are carried out under different con fining pressures(10 kPa,5 kPa and 2 kPa).Compared to DEM simulation,the simulated performance of SPH signi ficantly depends on the selected soil model and parameters.Only the appropriate soil model with optimal parameters adopted in SPH can generate the ideal simulation performance.To obtain a set of accurate parameters,the model parameters are determined through optimization-based or Bayesian-based parameter identi fication procedure(Jin et al.,2017b,2018b,2019b;Yin et al.,2017,2018b;Jin and Yin,2020).The parameters with more importance are put into the optimization procedure,as summarized in Table 6.An interesting phenomenon can be observed for the parameters optimized with the optimization-based or Bayesian-based parameter identi fication procedure.The largest difference in the model parameters is in the friction angle values,while other model parameters remain basically the same.One possible reason is that the different geometries(SS,TS and ES)lead to different contact conditions between the particles in the DEM simulation.As the particles become rough(from spherical to tetrahedral particles),the friction and shearing between the particles will increase on a microscopic level,resulting in the increase of friction angle on macroscopic level.Fig.4 presents the results of the SIMSAND model compared with the DEM triaxial data.A good agreement is observed between both simulations.

5.Results and discussion of slope collapse simulation

During the landslide process,the top width of the slope gradually decreases and the slip particles move along the upper boundary to the foot of the slope and slide to the far end.The maximum runout particle distance is namedDmand the final top width of the slopeLf.The“loose”slope generated by spherical particles that do not consider the particle rotation effect should have the largest landslide distance.In Fig.5,it can be seen that even in this extreme case,the maximum displacement of the FLAC and Lagrangian domains is very small(less than 1 mm),and it is therefore ignored.Only the large deformation area results,discretized with DEM or SPH,are discussed hereafter.

Figs.6 and 7 show the final failure modes of the 8 cases using the DEM and the SPH approaches,respectively.It can be seen from Fig.6 that the particle shape has a signi ficant effect on landslide morphology.The cases shown in Fig.6a and e have the largest runout distanceDm.The top of the slope is completely destroyed in the“loose”slope case(Fig.6e).For the cases of regular tetrahedron(Fig.6b)and elongated particles(Fig.6c),the slope runout distances are relatively small and the final deposit con figurations are similar,indicating that the mechanical properties of the two particle shapes are relatively close.For the cases shown in Fig.6d and h,where the anti-rotation effect of the particles is considered,the runout distanceDmdecreases and the top slope widthLfincreases,compared with the cases presented in Fig.6a and e.The results obtained by the SPH method(Fig.7)show the same trends.

In order to quantify the accuracy of the simulations,the final deposit con figurations provided by the DEMand SPH are presented in Fig.8.A good agreement can be observed for different soil densities and particle shapes.The values of the runout distanceDmand the slope widthLfby SPH follow closely the DEM benchmark results.

Fig.9 monitors the velocity field at several time steps for the 8 slopes considering different initial conditions provided by SPH.The runout distanceDmtends to be stable aroundT=3 s for all slopes.It can be observed that the particle velocity at the foot of all slopes is zero at this time step,only a small number of particles near the top of the slope are still moving.It is also obvious that the particle shape has a signi ficant effect on the landslide morphology.From Fig.9a-g and for different particle shapes(i.e.SS,TS,ES),the time required for the slope collapse to reach the final stable state(i.e.when the velocity of all particles reaches zero)decreases in sequence.Furthermore,as shown in Fig.9d and g,the rolling resistance increases the time for the collapse to reach the final stable state(compared with the SS cases where the rolling resistance is considered equal to zero).

Fig.5.Displacement field of a“loose”slope with spherical particles without rolling resistance:(a)DEM and(b)SPH.

Fig.6.Final failure patterns by DEM:(a)-(d)“Dense”slopes with(a)SS(no rolling resistance),(b)TS,(c)ES,and(d)SS(rolling resistance of 0.1);(e)-(h)“Loose”slopes with(e)SS(no rolling resistance),(f)TS,(g)ES,and(h)SS(rolling resistance of 0.1).

Fig.7.Final failure patterns by SPH:(a)-(d)“Dense”slopes with(a)SS(no rolling resistance),(b)TS,(c)ES,and(d)SS(rolling resistance of 0.1);(e)-(h)“Loose”slopes with(e)SS(no rolling resistance),(f)TS,(g)ES,and(h)SS(rolling resistance of 0.1).

Fig.8.Comparison of final deposit con figurations using DEMand SPH:(a)-(d)“Dense”slopes with(a)SS(no rolling resistance),(b)TS,(c)ES,and(d)SS(rolling resistance of 0.1);(e)-(h)“Loose”slopes with(e)SS(no rolling resistance),(f)TS,(g)ES,and(h)SS(rolling resistance of 0.1).

The porosity distributions of slope collapse at the final stable state are finally provided in Fig.10(SPH simulations).It can be observed that the porosity in the non-sliding area remains basically unchanged(before and after the landslide),while it increases in the sliding area.This is due to the particle redistribution after the landslide.For the loose slopes,it is also worth noting that the porosity distribution of the sliding area does not tend to be looser and the porosity does not signi ficantly change.One possible reason is that the spatial distribution of the particles after collapse is similar to that at the initial state.Another explanation could be that the shear expansion caused by the particles in the sliding area is not signi ficant.

6.Conclusions

Based on the critical state sand model SIMSAND and the SPH method,a highly ef ficient and accurate approach for numerical simulation of slope collapse is proposed in this paper.Using DEM results as benchmark studies,a total of 8 slopes with different initial conditions were simulated:3 different particle shapes(SS,TS and ES),and 2 different initial densities(dense and loose),and a set of numerical simulations considering the anti-rotation effect of spherical particles were studied.The results of landslide distance,particle rotation and velocity distribution were analyzed.

Comparisons show that the adopted numerical strategy is able to reproduce qualitatively and quantitatively the main characteristics of slope collapse,i.e.free surface,final deposit con figurations and final runout distance.The results also illustrate that the particle shape,initial density and particle anti-rotation coef ficient have a signi ficant effect on the landslide.Compared with regular spherical particles,the slope generated by irregular particles is less prone to sliding damage,its sliding area is smaller and the impact of landslide damage is less.Particles with irregular shape can therefore restrict the rotation of the particles.Compared with the“loose”slope,the“dense”slope presents less landslide damage.The numerical simulation results of regular tetrahedral and elongated particles are similar,indicating that these two particle shapes have similar restraining effects on particle rotation.

Fig.9.Velocity field variation during collapse using SPH:(a)-(d)“Dense”slopes with(a)SS(no rolling resistance),(b)TS,(c)ES,and(d)SS(rolling resistance of 0.1);(e)-(h)“Loose”slopes with(e)SS(no rolling resistance),(f)TS,(g)ES,and(h)SS(rolling resistance of 0.1).

The combination of SPH with the SIMSAND critical state constitutive model is able to reproduce slope collapse considering the effect of particle shapes and soil densities.The proposed strategy has a higher computational ef ficiency compared to DEM,while maintaining the computational accuracy.Therefore,it is an effective numerical tool that can be further applied to study real scale field landslides.

Fig.10.Porosity distribution at the final state using SPH:(a)-(d)“Dense”slopes with(a)SS(no rolling resistance),(b)TS,(c)ES,and(d)SS(rolling resistance of 0.1);(e)-(h)“Loose”slopes with(e)SS(no rolling resistance),(f)TS,(g)ES,and(h)SS(rolling resistance of 0.1).

Declarationofcompetinginterest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to in fluence the work reported in this paper.

Acknowledgments

This research is supported by Shenzhen(China)Science and Technology Innovation Committee (Grant Nos.JSGG2018 0504170449754),and the numerical simulations are supported by Center for Computational Science and Engineering at Southern University of Science and Technology,Shenzhen,China.