Novel SPH SIMSAND–Based Approach for Modeling of Granular Collapse

Granular collapse is a common issue in natural hazards. This paper proposes a novel numerical approach on modeling granular column collapse. A newly developed critical state-based constitutive model, SIMSAND, was adopted to combine with the smoothed particle hydrodynamics (SPH) method for realistically reproducing large deformation during collapse. A rectangular channel and two-dimensional column tests were first simulated for the validation. The effects of aspect ratio and initial soil density were further investigated by additional simulations. It was demonstrated that the novel SPH-SIMSAND approach is helpful in improving the understanding of granular collapse and should be an effective computational tool for the analysis of real-scale granular flow.


Introduction
Granular collapse, such as debris, rock avalanches, and landslides, is a common issue in natural hazards. To understand this phenomenon, two types of experiments have been performed: rectangular channel flowing tests (Balmforth and Kerswell 2005;Lajeunesse et al. 2005;Lube et al. 2005Lube et al. , 2007Bui et al. 2008;Crosta et al. 2009) and column flowing tests (Lube et al. 2004(Lube et al. , 2005(Lube et al. , 2007. In the first, granular collapse is obtained by putting the granular material in a rectangular channel and quickly removing a vertical side boundary, whereas in the second, the granular material is in a hollow cylinder pipe. Daerr and Douady (1999) demonstrated that the final deposit morphology for the column collapse tests depended on the initial soil density based on only a low aspect ratio test. Later, experimental results by others (Lajeunesse et al. 2004;Lube et al. 2004Lube et al. , 2005Lube et al. , 2007 showed that the final deposit morphology (deposit radius, deposit height, and slumping velocity) mainly depended on the initial aspect ratio of the granular column without showing the influence of soil density. More recently, various numerical studies were also conducted to investigate the granular collapse using the physics-based discrete-element method (DEM). The results also showed that the final deposit morphology depended on the initial aspect ratio (Staron and Hinch 2005;Zenit 2005;Lacaze et al. 2008;Girolami et al. 2012;Soundararajan 2015;Utili et al. 2015). In these studies, the effect of the initial void ratio (corresponding to the soil density) was again ignored. To confirm this, Kermani et al. (2015) and Soundararajan (2015) simulated the effect of initial porosity on three-dimensional (3D) asymmetrical collapse considering different aspect ratios. However, the number of particles in the majority of DEM simulations is limited and far from a real physical model or case due to computational efficiency issues. The applicability of the DEM to a real-scale problem is, thus, still questionable.
The finite-element method is a powerful way to simulate problems in geotechnics (Shen et al. 2014Wu et al. 2016Wu et al. , 2017a. The method has also been applied to analyze granular column collapse using Mohr-Coulomb and Drucker-Prager constitutive models and the arbitrary Lagrangian-Eulerian (ALE) technique (Crosta et al. 2009), the particle finite-element method (PFEM) (Zhang et al. 2015), the smooth particle hydrodynamics (SPH) method (Bui et al. 2008;Nguyen et al. 2017), and the material point method (MPM) (Sołowski and Sloan 2015). However, these constitutive models adopted up to now are not appropriate for describing the state-dependent behavior of soils. Because the constitutive model is a key component controlling the physics of finite-element analysis, the well-accepted critical-state modeling theory for granular material should be considered.
In this paper, a critical state-based model accounting for soil density effects was first adopted and implemented in combination with the SPH method for large-deformation analysis. Then, rectangular channel soil collapse tests were simulated and the combined numerical tool was validated. Granular column collapse tests with different aspect ratios were simulated for further validation in terms of the final deposit morphologies, dynamic flowing profiles, and undisturbed areas. Finally, additional simulations for different initial void ratios (e 0 = 0.75, 0.85, 0.95, and 1.05) were conducted to investigate the effect of soil density on the final deposit morphology. The collapse evolution process was monitored, and the ability of the adopted modeling strategy to reproduce the influences of the initial aspect ratio and the soil density was evaluated.

SPH Method
The SPH method implemented in the finite-element code ABAQUS/ Explicit was adopted in this study to solve large-deformation problems. The method was initially proposed by Gingold and Monaghan (1977) for numerical analysis in astrophysics. Further developments allowed for applications to various problems in the field of solid mechanics. In the SPH method, number of particles is used to discretize the numerical domain. The particle represents the volume and mass of soils and carries numerical parameters, such as acceleration, velocity, and void ratio (Chen and Qiu 2012), as shown in Fig. 1(a).
The field variable fx ðÞat point x in the numerical domain can then be obtained by the following equation accounting for the effect of neighboring particles: where W = weighting function, or the so-called kernel or smoothing function. The fx ðÞis further approximated by the summation over neighboring particles by fx ðÞ¼ where V i , r i , and m i = volume, density, and mass of the particle i, respectively; and N = number of influencing particles. The spatial derivative of fx ðÞis approximated by the differential operations on the kernel function Thus, the kernel function has a significant influence on the efficiency and accuracy of SPH calculations. The particles serve as interpolation points for estimating all variables in the continuous medium. These particles can also be separated by a big distance. The variables between particles can be approximated by smoothing shape functions. When a particle gets to a certain distance, which is the smoothing length (h), from another one, particles start to interact. Only when the SPH particles are in the influence domain will they interact with each other. Thus, a smoother or more continuous behavior can be achieved by a larger smoothing length. Otherwise, more discrete behaviors will be obtained using a smaller smoothing length, because particles are more independent in this case.
The main advantage of the SPH method is that, when calculating spatial derivatives, there is no need for a fixed computational grid. Instead, using analytical expressions based on the derivatives of the smoothing functions, estimates of derivatives canbeobtained (Li and Liu 2002). Because the particles during the calculation can be interacted and separated over the time, the SPH method can deal with very-large-deformation analysis (Bojanowski 2014).

Explicit Finite-Element Method
The SPH method implemented in ABAQUS uses the explicit timeintegration method (Hibbitt, Karlsson and Sorensen 2001). As shown in Fig. 1(b), the equilibrium condition is first written with the balance of internal force and external force where M = mass matrix; € u = acceleration; P = applied external force vector; and I = internal force vector. The equations of motion for the body are then integrated using the explicit time central-difference integration rule as follows: where u = displacement; subscript t = time in an explicit dynamic step; _ u = velocity; and Dt = increment of time. For the stability of calculation, the time increment (Dt) should be smaller than a limited value (Dt % L min =c d ) with the smallest element dimension of mesh (L min ) and the dilatory wave speed (c d ).
Finally, incremental strain calculated by incremental displacement is used by the constitutive model to update stresses and then internal forces, up to a new equilibrium condition.

Adopted Critical State-Based Model SIMSAND
The SIMSAND model was developed based on the Mohr-Coulomb model by implementing the critical-state concept (Yin et al. 2017a;Jin et al. 2017) with nonlinear elasticity, nonlinear plastic hardening, and a simplified 3D strength criterion. The state-dependent peak strength and stress dilatancy (contraction or dilation) are well captured by the SIMSAND model (Jin et al. 2017). The basic constitutive equations are summarized in Table 1. The model parameters with their definitions are summarized in Table 2. The calibration of the model parameters can be carried out in a straightforward way (Wu et al. 2017b) or using optimization methods (Jin et al. 2016a, b;Yin et al. 2017b). The adopted SIMSAND model was implemented into ABAQUS/ Explicit as a user-defined material model via user material subroutine VUMAT. The procedure of model implementation followed that of Hibbitt, Karlsson and Sorensen (2001). In ABAQUS/Explicit combined with VUMAT, the strain increment on the element D« at Dt was first solved by ABAQUS using the presented explicit time central-differential integration method. Then, the stress increment (Ds ) was updated through VUMAT using the solved D« . For the stress integration, the cutting-plane algorithm proposed by Ortiz and Simo (1986)was adopted. Implementation of the finite-element method was then verified by simulating the aforementioned drained triaxial tests before modeling the granular collapse.

Experimental and Numerical Simulation Configurations
To validate the SIMSAND model and the adopted numerical integration scheme, the rectangular channel soil collapse tests (Bui et al. 2008) were simulated hereafter. In the experiments, small aluminum bars of various diameters (0.1 and 0.15 cm) were used to model the soil. The bars were initially arranged into an area of 20cm length Â 10-cm height Â 2-cm width, delimited by two flat solid walls. The experiment was started by quickly moving the right wall horizontally to the right, causing the flow of the aluminum bars to the side due to gravity (Fig. 2).
For the numerical model, the spatial discretization domain is shown in Fig. 3(a). SPH particles were used to model the soil while the two solid walls were discretized with rigid hexahedral finite elements. The initial SPH particle distance was (approximately) the same in the horizontal and vertical directions to reproduce homogeneous conditions.
A cell size of 0.2 cm was estimated by checking mesh dependency (Fig. 3 with different size from 0.1 to 0.2 cm) so that the particle space or the particle diameter was half of the cell size (eight particles in each cell). The total number of particles was around 400,000. The bottom surface was set as a fixed boundary, whereas symmetrical conditions were assumed for the four lateral boundaries.
The simulation was carried out in two steps: to balance the geostatic field and to move the right wall along the horizontal direction with a speed of 1 m/s to the right. The contacts between the sand Critical state with interlocking effects sin 3u "# 1 4 with c 2 ¼ 3 À sin f pt 3 þ sin f pt Note: p at = atmospheric pressure (p at = 101.3 kPa); p 0 = effective mean pressure; q = deviatoric stress; e c = critical-state void ratio; f p = peak friction angle; f pt = phase transformation friction angle; M pt = stress ratio corresponding to the phase transformation; M p = peak stress ratio; u = Lode's angle with its effect introduced by using Sheng et al. (2000), which is similar to Yao et al. (2004Yao et al. ( , 2008Yao et al. ( , 2009; and n p and n d = interlocking parameters controlling the degree of interlocking due to neighboring particles according to Yin et al. (2010) and Yin and Chang (2013). Note: CSL = Critical State Line. and walls were described by the classical Coulomb friction law with a friction coefficient ¼ tan(f /2) ¼ 0.28. Following Lube et al. (2004Lube et al. ( , 2005, the collapse behavior was less sensitive to the grain properties of granular material compared to the density. Meanwhile, due to a lack of mechanical tests on used material, the model parameters of Toyoura sand were adopted for simulations. Fig. 4 presents the typical statedependent behavior of granular materials simulated by the model under a very low confining stress in the experiments. As showninFig.4, the initial void ratio (or initial density) significantly influenced the peak strength and the dilation/contraction of granular material, which should be considered in granular collapse because the density changes greatly during collapse.

Numerical Validation
Different void ratios from 0.75 to 1.05 were assumed for simulations and compared with experimental results. As shown in Fig. 5, the higher the initial density (or the smaller the void ratio), the stiffer the free surface and failure line were and the longer the runout distance was. This agrees with the experimental results of Lube et al. (2005). Luckily, the simulation with e 0 =0 . 9 5fit the experiment well. Thus, e 0 = 0.95 was adopted as a reference void ratio for the following sections.
During simulation using the SPH method, each particle represents one gauss integration point. Accordingly, similar to the element in FEM, the total strain of each particle is divided into the elastic and plastic strains when using the SIMSAND model in SPH. In this study, the equivalent plastic strain was defined as ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2=3 _ e p ij : _ e p ij r (where _ e p ij is the tensor of the deviatoric plastic strain rate) and was used to describe the plastic deformation. Fig. 6 shows the calculated deformed column with the distribution of deviatoric plastic strain at different time steps. The movement took place for approximately 0.6 s.

Numerical Simulation Configurations
To have a further understanding of the granular column collapse, the two-dimensional column collapse with a different initial aspect ratio was simulated.
In the numerical model, the spatial discretization domain had the same dimensions as in the experiments by Lube et al. (2005), as shown in Fig. 7, where h i is the initial height, d i is the initial basal length, and a = h i /d i is the aspect ratio of the granular column. As in the experiments, six aspect ratios were studied (a = 0.5, 1.0, 1.5, 3.0, 7.0, and 9.0) in which the column initial basal length was constant and equal to 10 cm. The element size for all aspect ratios was 0.2 cm, with which a total number of elements from 31,250 to 562,500 were discretized for different aspect ratios of columns, as summarized in Table 3.

Comparisons and Discussions
Three different aspect ratios (a = 0.5, 1.5, and 7.0) in the numerical results were compared with the experimental results of Lube et al. (2005)inFi g.8, where the gradually varied shading represents the distribution of the deviatoric plastic strain. A good agreement was observed in terms of deposit morphology. Furthermore, simulations captured the progressive collapsing process; more specifically: 1. For a = 0.5, the outer region at the bottom flowed and the length of the runout region at the column foot increased. The inner region was less disturbed than other cases, and a flat area remained at the top. 2. For a = 1.5, the runout distance at the bottom increased and a flat undisturbed area was created at the top during the initial stage. A cone tip formed at the end. 3. For a = 7.0, an important degradation of the column height occurred first. Then, the runout distance at the column foot increased and the upper initially undisturbed surface began to flow. Finally, a very large runout distance with a cone tip formed.  Flow Description Fig. 9 shows the successive numerical granular collapse profiles for the three aspect ratios. Lube et al. (2004) summarized three deposit morphologies based on an aspect ratio range: (1) a < 0.74, (2) 0.74 < a < 1.7, and (3) a > 1.7. These distinctive flow processes were well captured by simulations: 1. For a = 0.5 (<0.74), a lateral flow developed at the column foot, and a flat undisturbed area remained at the top. The deposit height (h 0 ) stayed constant. 2. For 0.74 < a = 1.5 < 1.7, the evolution of the lateral flow was accompanied by a small decrease of the initial height (h 0 ). A wedge shape formed at the end. 3. For a =7> 1.7, initially, the column height greatly decreased but the upper surface remained unchanged. Then, the lateral flow developed quickly. Simultaneously, the length of the upper surface decreased and formed a dome-like shape. At the final stage, the runout distance (d 1 Þ was important, and a wedge shape formed at the top h 1 ð ). Fig. 10 shows the comparison of plastic strain fields between simulations with different aspect ratios, where the smallest value of deviatoric plastic strain (ɛ p d ) is black, and the undisturbed stable area inside the granular column is suggested by a relatively small value of deviatoric plastic strain. Then, it was found that an undisturbed trapezoid area developed on the upper free surface of the column only for the cases of small aspect ratios (a = 0.5, 1.0, and 1.5), highlighted in black in Fig. 10. However, for the cases of large aspect ratios (a = 3.0, 7.0, and 9.0), the upper free surface developed a triangle area and the repose angle presented an increasing trend.

Deposit Morphology
To study the influence of soil density, granular columns with four initial void ratios (e 0 = 0.75, 0.85, 0.95, and 1.05, corresponding to unit weights of g = 1.51, 1.43, 1.36, and 1.29 when G s = 2.65) and six aspect ratios (a = 0.5, 1.0, 1.5, 3.0, 7.0, and 9.0) were simulated hereafter. The numerical results were compared with the best-fitting equations shown in Fig. 11. The numerical results were in agreement with the best-fitting equations of Lube et al. (2005) for the normalized final runout distance when e 0 = 0.95. However, a difference appeared for smaller initial void ratios, especially in the cases of larger aspect ratios. The normalized final deposit height seemed less sensitive to the soil density. However, the effect increased for the cases of higher aspect ratios. All comparisons indicated that the deposit morphology (final runout distance and final deposit height) depended on not only the aspect ratio but also the initial density. This was also in agreement with DEM results (Kermani et al. 2015). Fig. 12 presents the final deposit morphology of the simulations with different initial void ratios. The final deposit morphology was sensitive to the initial void ratio or soil density. For the cases of lower column a < 1.15, the height of the circular undisturbed zone in the upper surface remained the same, but its area decreased. The runout distance increased with the increase of the initial void ratio (or decrease of initial soil density). For the cases of higher column a > 1.15, the denser granular material caused a shorter runout distance and higher final deposited height in agreement with the DEM simulations by Kermani et al. (2015) and experiments by Daerr and Douady (1999). This illustrates that the influence of the void ratio or soil density on granular collapse was well captured by using the SPH technique with the critical state-based SIMSAND model. The difference in deposit morphology between different void ratios for the same aspect ratio can be explained by the stress dilatancy of granular soils. For dense sand, a stronger interlocking force between particles developed with a higher mobilized strength and finally formed a stable inner region during the stage of collapsing. Therefore, a denser granular column corresponded to a bigger inclination of the slope surface with a higher deposit height and smaller runout distance.

Monitoring Collapse
The collapsing time (t) was normalized by the intrinsic critical time (t c ) used for the dimensionless analysis, where the value of t c can be calculated by the initial height of the sand column (t c = ffiffiffiffiffiffiffiffi ffi h i =g p ) (Soundararajan 2015   results present an S-shaped curve with two successive stages, regardless of the initial aspect ratios and void ratios. First there was an acceleration and then a deceleration stage starting close to 1.5t c -2.0t c . The collapse ceased at approximately 3.5t c . The following conclusions can also be made: 1. For the same initial void ratio, a higher aspect ratio leads to a more important normalized runout distance and deposit height. 2. For the same aspect ratio, denser sand (smaller initial void ratio) leads to a more important deposit height but to a shorter runout distance in the same timescale. 3. The normalized deposit height initially increases when the right boundary plate is lifted [ Figs. 14(a-c)]. This uplift depends strongly on the soil density. A denser soil leads to a more important uplift. The reason for this is again the density effect, which develops more on denser sand due to the stronger interlocking between particles at initial shearing.

Conclusion
A numerical investigation on granular collapse based on the critical-state soil model SIMSAND and the SPH method was carried out. The validation was first provided by comparing experimental data from rectangular channel and two-dimensional column   collapse tests. Then, the influence of the initial aspect ratio and the soil density was studied in detail. All comparisons showed that the adopted numerical strategy was able to reproduce qualitatively and quantitatively the main behaviors of granular column collapse (i.e., free surface, failure line, final deformed profile for the rectangular channel test, final runout distance, and deposit height). More specifically, when the initial soil density decreased, the failure surface shrank and the free surface enlarged. A lower initial void ratio generated a stronger interlocking force, leading to a higher deposit height and a shorter runout distance.    The combination of the SIMSAND model with the SPH method is able to reproduce granular collapse, taking into account the influence of different aspect ratios and soil densities. Therefore, it provides an effective computational tool for the analysis of real-scale granular flow.