# Novel SPH SIMSAND-Based Approach for Modeling of Granular Collapse

Zhen-Yu Yin, Zhuang Jin, Panagiotis Kotronis, Ze-Xiang Wu

► To cite this version:

Zhen-Yu Yin, Zhuang Jin, Panagiotis Kotronis, Ze-Xiang Wu. Novel SPH SIMSAND-Based Approach for Modeling of Granular Collapse. International Journal of Geomechanics, 2018, 18 (11), 10.1061/(ASCE)GM.1943-5622.0001255. hal-01876908

HAL Id: hal-01876908

https://hal.science/hal-01876908

Submitted on 24 Jun 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Novel SPH SIMSAND-Based Approach for Modeling of
Granular Collapse

Zhen-Yu Yin<fnref n="1" />; Zhuang Jin<fnref n="2" />; Panagiotis Kotronis<fnref n="3" />; and Ze-Xiang Wu<fnref n="4" />

**Abstract:** 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.

**keywords:** Granular material; Collapse; Smoothed particle hydrodynamics (SPH); Large-deformation analysis; Granular flow; Critical state.

## 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. 2005, 2007; Bui et al. 2008; Crosta et al. 2009) and column flowing tests (Lube et al. 2004, 2005, 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. 2004, 2005, 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; Utiли 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. 2014, 2017; Wu et al. 2016, 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.

<footnote n="1">¹Associate Professor, Dept. of Civil and Environmental Engineering, Hong Kong Polytechnic Univ., Hung Hom, Kowloon, Hong Kong, China (corresponding author). Email: zhenyu.yin@gmail.com; zhenyu.yin@polyu.edu.hk</footnote>

<footnote n="2">²Ph.D. Candidate, Research Institute of Civil Engineering and Mechanics (GeM), UMR CNRS 6183, Ecole Centrale de Nantes, Nantes 44321, France.</footnote>

<footnote n="3">³Professor, Research Institute of Civil Engineering and Mechanics (GeM), UMR CNRS 6183, Ecole Centrale de Nantes, Nantes 44321, France.</footnote>

<footnote n="4">⁴Assistant Professor, College of Architecture and Civil Engineering, Wenzhou Univ., Wenzhou 325035, People's Republic of China.</footnote>
# Numerical Methods and Constitutive Model for Simulation

## 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 $f(x)$ at point $x$ in the numerical domain can then be obtained by the following equation accounting for the effect of neighboring particles:

$$f(x) = \int_{\Omega} f(x')W(x - x', h)dx' \quad (1)$$

where $W$ = weighting function, or the so-called kernel or smoothing function. The $f(x)$ is further approximated by the summation over neighboring particles by

$$f(x) = \sum_{i=1}^{N} f(x_i)W(x - x_i, h)V_i = \sum_{i=1}^{N} f(x_i)W(x - x_i, h) \frac{m_i}{\rho_i} \quad (2)$$

where $V_i$, $\rho_i$, and $m_i$ = volume, density, and mass of the particle $i$, respectively; and $N$ = number of influencing particles. The spatial derivative of $f(x)$ is approximated by the differential operations on the kernel function

$$\frac{\partial f(x)}{\partial x} = \sum_{i=1}^{N} \frac{m_i}{\rho_i} f(x_i) \frac{\partial W(x - x_i, h)}{\partial x_i} \quad (3)$$

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 can be obtained (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 time-integration 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

$$M\ddot{u}_{(t)} = P_{(t)} - I_{(t)} \quad (4)$$

where $M$ = mass matrix; $\ddot{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:

$$\begin{aligned} \dot{u} \left( t + \frac{\Delta t}{2} \right) &= \dot{u} \left( t - \frac{\Delta t}{2} \right) + \frac{\Delta t_{(t+\Delta t)} + \Delta t_{(t)}}{2} \ddot{u}_{(t)} \\ u_{(t+\Delta t)} &= u_{(t)} + \Delta t_{(t+\Delta t)} \dot{u} \left( t + \frac{\Delta t}{2} \right) \end{aligned} \quad (5)$$

where $u$ = displacement; subscript $t$ = time in an explicit dynamic step; $\dot{u}$ = velocity; and $\Delta t$ = increment of time. For the stability of calculation, the time increment ($\Delta t$) should be smaller than a limited value ($\Delta t \approx 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.

<figure><img src="image_2.png" /><figcaption>Fig. 1. SPH numerical method: (a) deformation of SPH method; and (b) flowchart of explicit analysis.</figcaption></figure>
## 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 sub-routine 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 $\Delta\varepsilon$ at $\Delta t$ was first solved by ABAQUS using the presented explicit time central-differential integration method. Then, the stress increment ($\Delta\sigma$) was updated through VUMAT using the solved $\Delta\varepsilon$. 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.

## Rectangular Channel Soil Collapse Simulations

### 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 20-cm 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

<table><caption>**Table 2.** Parameters of SIMSAND model and values of parameters for Toyoura sand</caption><thead><tr><th>Parameters</th><th>Definitions</th><th>Value</th></tr></thead><tbody><tr><td>K<sub>0</sub></td><td>Referential bulk modulus (dimensionless)</td><td>130</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>0.2</td></tr><tr><td>n</td><td>Elastic constant controlling nonlinear stiffness</td><td>0.5</td></tr><tr><td>φ</td><td>Critical-state friction angle</td><td>31.5</td></tr><tr><td>e<sub>ref</sub></td><td>Initial critical-state void ratio</td><td>0.977</td></tr><tr><td>λ</td><td>Constantly controlling the nonlinearity of CSL</td><td>0.0596</td></tr><tr><td>ξ</td><td>Constantly controlling the nonlinearity of CSL</td><td>0.365</td></tr><tr><td>A<sub>d</sub></td><td>Constant of magnitude of the stress dilatancy</td><td>0.7</td></tr><tr><td>k<sub>p</sub></td><td>Plastic modulus-related constant</td><td>0.0044</td></tr><tr><td>n<sub>p</sub></td><td>Interlocking related peak strength parameter</td><td>2.4</td></tr><tr><td>n<sub>d</sub></td><td>Interlocking related phase transformation parameter</td><td>2.9</td></tr></tbody></table>

Note: CSL = Critical State Line.

<table><caption>**Table 1.** Basic constitutive equations of SIMSAND model</caption><thead><tr><th>Components</th><th>Constitutive equations</th></tr></thead><tbody><tr><td>Elasticity</td><td>ε̇<sub>ij</sub><sup>e</sup> = (1 + ν) σ' <sub>ij</sub> - (ν / 3K(1 - 2ν)) σ' <sub>kk</sub> δ<sub>ij</sub> K = K<sub>0</sub> · p<sub>at</sub> (2.97 - e)<sup>2</sup> (p' / p<sub>at</sub>)<sup>ξ</sup></td></tr><tr><td>Yield surface</td><td>f = q / p' - H</td></tr><tr><td>Potential surface</td><td>∂g / ∂p' = A<sub>d</sub> (M<sub>pt</sub> - q / p'); ∂g / ∂s<sub>ij</sub> = { 1 1 1 1 1 1 1 }</td></tr><tr><td>Hardening rule</td><td>H = M<sub>p</sub>ε<sub>d</sub><sup>p</sup> / (k<sub>p</sub> + ε<sub>d</sub><sup>p</sup>)</td></tr><tr><td>Critical state with interlocking effects</td><td>e<sub>c</sub> = e<sub>ref</sub> - λ (p' / p<sub>at</sub>)<sup>ξ</sup> tan φ<sub>p</sub> = (e<sub>c</sub> / e)<sup>n<sub>p</sub></sup> tan φ<sub>μ</sub> tan φ<sub>pt</sub> = (e<sub>c</sub> / e)<sup>-n<sub>d</sub></sup> tan φ<sub>μ</sub></td></tr><tr><td>3D strength criterion</td><td>M<sub>p</sub> = (6 sin φ<sub>p</sub> / (3 - sin φ<sub>p</sub>)) [2c<sub>1</sub><sup>4</sup> / (1 + c<sub>1</sub><sup>4</sup> + (1 - c<sub>1</sub><sup>4</sup>)sin 3θ)]<sup>1/4</sup> with c<sub>1</sub> = (3 - sin φ<sub>p</sub> / (3 + sin φ<sub>p</sub>))<br>M<sub>pt</sub> = (6 sin φ<sub>pt</sub> / (3 - sin φ<sub>pt</sub>)) [2c<sub>2</sub><sup>4</sup> / (1 + c<sub>2</sub><sup>4</sup> + (1 - c<sub>2</sub><sup>4</sup>)sin 3θ)]<sup>1/4</sup> with c<sub>2</sub> = (3 - sin φ<sub>pt</sub> / (3 + sin φ<sub>pt</sub>))</td></tr></tbody></table>

Note: $p_{at}$ = atmospheric pressure ($p_{at}$ = 101.3 kPa); $p'$ = effective mean pressure; $q$ = deviatoric stress; $e_c$ = critical-state void ratio; $\phi_p$ = peak friction angle; $\phi_{pt}$ = phase transformation friction angle; $M_{pt}$ = stress ratio corresponding to the phase transformation; $M_p$ = peak stress ratio; $\theta$ = Lode's angle with its effect introduced by using Sheng et al. (2000), which is similar to Yao et al. (2004, 2008, 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).
and walls were described by the classical Coulomb friction law with a friction coefficient = tan($\phi/2$) = 0.28.

Following Lube et al. (2004, 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 state-dependent behavior of granular materials simulated by the model under a very low confining stress in the experiments. As shown in Fig. 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.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Rectangular channel soil collapse experiments.</figcaption></figure>

## 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.95 fit 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 $\sqrt{2/3(\dot{e}_{ij}^p : \dot{e}_{ij}^p)}$ (where $\dot{e}_{ij}^p$ 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.

## Influence of Aspect Ratio

### 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,

<figure><img src="image_4.png" /><figcaption>Fig. 3. SPH particle density and mesh dependency: (a) two-dimensional (2D) spatial domain discretization for SPH particles and rigid hexahedral finite elements; (b) element size of 0.1 cm; (c) element size of 0.15 cm; and (d) element size of 0.2 cm.</figcaption></figure>
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) in Fig. 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.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Simulated stress–strain behavior of conventional triaxial tests using parameters of Toyoura sand: (a) deviatoric stress versus axial strain; (b) void ratio versus axial strain; and (c) void ratio versus mean effective stress in demilogarithmic scale. CSL = Critical State Line.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Simulation results of column collapse: (a) final free surfaces and failure lines for various initial void ratios $e_0$ = 0.75, 0.85, 0.95, 1.05; (b) final deformed column in experiment; and (c) final deformed column by simulation with $e_0$ = 0.95.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Deformed shape and deviatoric plastic strain distribution of the granular column at different time steps.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Discretization domain plane views for six aspect ratios (a).</figcaption></figure>

<table><caption>**Table 3.** Discretization parameters for all simulations</caption><thead><tr><th>Aspect ratio (a)</th><th>Initial basal length [d<sub>i</sub> (cm)]</th><th>Initial height [h<sub>i</sub> (cm)]</th><th>Column size (cm<sup>3</sup>)</th><th>Number of SPH particles</th></tr></thead><tbody><tr><td>0.5</td><td>10</td><td>5</td><td>10 × 5 × 5</td><td>31,250</td></tr><tr><td>1.0</td><td>10</td><td>10</td><td>10 × 10 × 5</td><td>62,500</td></tr><tr><td>1.5</td><td>10</td><td>15</td><td>10 × 15 × 5</td><td>93,750</td></tr><tr><td>3.0</td><td>10</td><td>30</td><td>10 × 30 × 5</td><td>187,500</td></tr><tr><td>7.0</td><td>10</td><td>70</td><td>10 × 70 × 5</td><td>437,500</td></tr><tr><td>9.0</td><td>10</td><td>90</td><td>10 × 90 × 5</td><td>562,500</td></tr></tbody></table>

## 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_\infty$) was important, and a wedge shape formed at the top ($h_\infty$).

Fig. 10 shows the comparison of plastic strain fields between simulations with different aspect ratios, where the smallest value of deviatoric plastic strain ($\epsilon_d^p$) 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.

## Influence of Soil Density

### 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 $\gamma = 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 = \sqrt{h_i/g}$) (Soundararajan 2015). Figs. 13 and 14 present the simulation results of the normalized runout distance $[(d_\infty - d_i)/d_i]$ and the normalized deposit height ($h_\infty/d_i$) with the normalized time ($t/t_c$). All
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison of progressive granular collapse between simulations and experiments for different aspect ratios of (a) *a* = 0.5; (b) *a* = 1.5; and (c) *a* = 7.0.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Numerical progressive column collapse profiles for different aspect ratios of (a) *a* = 0.5; (b) *a* = 1.5; and (c) *a* = 7.0.</figcaption></figure>

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
<figure><img src="image_11.png" /><figcaption>Fig. 10. Side view of plastic deviatoric strain after collapse with different a: (a) a = 0.5; (b) a = 1.0; (c) a = 1.5; (d) a = 3.0; (e) a = 7.0; and (f) a = 9.0.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Comparison between the numerical simulations and the best-fitting equations of Lube et al. (2005) for different aspect ratios and initial void ratios: (a) normalized final runout distance; and (b) normalized final deposit height.</figcaption></figure>

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.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 12.** Influence of the initial void ratio on the final deposit morphology for (a) $a = 0.5$; (b) $a = 1.0$; (c) $a = 1.5$; (d) $a = 3.0$; (e) $a = 7.0$; and (f) $a = 9.0$.</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 13.** Evolution of the normalized runout distance for different aspect ratios and initial void ratios: (a) for $a = 0.5$, (b) for $a = 1.5$; and (c) for $a = 7.0$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 14. Evolution of the normalized deposit height for different aspect ratios and initial void ratios: (a) for $a = 0.5$, (b) for $a = 1.5$; and (c) for $a = 7.0$.</figcaption></figure>

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.

## Acknowledgments

The financial support for this research came from the National Natural Science Foundation of China (Grants 41372285 and 51579179) and the Region Pays de la Loire of France (Project RI-ADAPTCLIM).

## References

* Balmforth, N., and R. Kerswell. 2005. "Granular collapse in two dimensions." *J. Fluid Mech.* 538 (1): 399-428. https://doi.org/10.1017/S0022112005005537.

* Bojanowski, C. 2014. "Numerical modeling of large deformations in soil structure interaction problems using FE, EFG, SPH, and MM-ALE formulations." *Arch. Appl. Mech.* 84 (5): 743-755. https://doi.org/10.1007/s00419-014-0830-5.

* Bui, H. H., R. Fukagawa, K. Sako, and S. Ohno. 2008. "Lagrangian mesh-free particles method (SPH) for large deformation and failure flows of geomaterial using elastic-plastic soil constitutive model." *Int. J. Numer. Anal. Methods Geomech.* 32 (12): 1537-1570. https://doi.org/10.1002/nag.688.

* Chen, W., and T. Qiu. 2012. "Numerical simulations for large deformation of granular materials using smoothed particle hydrodynamics method." *Int. J. Geomech.* 12 (2): 127-135. https://doi.org/10.1061/(ASCE)GM.1943-5622.0000149.

* Crosta, G., S. Imposimato, and D. Roddeman. 2009. "Numerical modeling of 2-D granular step collapse on erodible and nonerodible surface." *J. Geophys. Res. Earth Surf.* 114 (F3): F03020. https://doi.org/10.1029/2008JF001186.

* Daerr, A., and S. Douady. 1999. "Sensitivity of granular surface flows to preparation." *Europhys. Lett.* 47 (3): 324-330. https://doi.org/10.1209/epl/1999-00392-7.

* Gingold, R. A., and J. J. Monaghan. 1977. "Smoothed particle hydrodynamics: Theory and application to non-spherical stars." *Mon. Not. R. Astron. Soc.* 181 (3): 375-389. https://doi.org/10.1093/mnras/181.3.375.

* Girolami, L., V. Hergault, G. Vinay, and A. Wachs. 2012. "A three-dimensional discrete-grain model for the simulation of dam-break reticular collapses: Comparison between numerical results and experiments." *Granular Matter* 14 (3): 381-392. https://doi.org/10.1007/s10035-012-0342-3.

* Hibbitt, Karlsson and Sorensen. 2001. *ABAQUS/explicit: User's manual*. Hibbitt, Karlsson and Sorensen Incorporated.

* Jin, Y.-F., Z.-X. Wu, Z.-Y. Yin, and J. S. Shen. 2017. "Estimation of critical state-related formula in advanced constitutive modeling of granular

material." *Acta Geotech.* 12 (6): 1329-1351. https://doi.org/10.1007/s11440-017-0586-5.

* Jin, Y.-F., Z.-Y. Yin, S.-L. Shen, and P.-Y. Hicher. 2016a. "Investigation into MOGA for identifying parameters of a critical-state-based sand model and parameters correlation by factor analysis." *Acta Geotech.* 11 (5): 1131-1145. https://doi.org/10.1007/s11440-015-0425-5.

* Jin, Y. F., Z. Y. Yin, S. L. Shen, and P. Y. Hicher. 2016b. "Selection of sand models and identification of parameters using an enhanced genetic algorithm." *Int. J. Numer. Anal. Methods Geomech.* 40 (8): 1219-1240. https://doi.org/10.1002/nag.2487.

* Kermani, E., T. Qiu, and T. Li. 2015. "Simulation of collapse of granular columns using the discrete element method." *Int. J. Geomech.* 15 (6): 04015004. https://doi.org/10.1061/(ASCE)GM.1943-5622.0000467.

* Lacaze, L., J. C. Phillips, and R. R. Kerswell. 2008. "Planar collapse of a granular column: Experiments and discrete element simulations." *Phys. Fluids* 20 (6): 063302. https://doi.org/10.1063/1.2929375.

* Lajeunesse, E., A. Mangeney-Castelnau, and J. Vilotte. 2004. "Spreading of a granular mass on a horizontal plane." *Phys. Fluids* 16 (7): 2371-2381. https://doi.org/10.1063/1.1736611.

* Lajeunesse, E., J. Monnier, and G. Homsy. 2005. "Granular slumping on a horizontal surface." *Phys. Fluids* 17 (10): 103302. https://doi.org/10.1063/1.2087687.

* Li, S., and W. K. Liu. 2002. "Meshfree and particle methods and their applications." *Appl. Mech. Rev.* 55 (1): 1-34. https://doi.org/10.1115/1.1431547.

* Lube, G., H. E. Huppert, R. S. J. Sparks, and A. Freundt. 2005. "Collapses of two-dimensional granular columns." *Phys. Rev. E* 72 (4): 041301. https://doi.org/10.1103/PhysRevE.72.041301.

* Lube, G., H. E. Huppert, R. S. J. Sparks, and A. Freundt. 2007. "Static and flowing regions in granular collapses down channels." *Phys. Fluids* 19 (4): 043301. https://doi.org/10.1063/1.2712431.

* Lube, G., H. E. Huppert, R. S. J. Sparks, and M. A. Hallworth. 2004. "Axisymmetric collapses of granular columns." *J. Fluid Mech.* 508: 175-199. https://doi.org/10.1017/S0022112004009036.

* Nguyen, C. T., C. T. Nguyen, H. H. Bui, G. D. Nguyen, and R. Fukagawa. 2017. "A new SPH-based approach to simulation of granular flows using viscous damping and stress regularisation." *Landslides* 14 (1): 69-81. https://doi.org/10.1007/s10346-016-0681-y.

* Ortiz, M., and J. Simo. 1986. "An analysis of a new class of integration algorithms for elastoplastic constitutive relations." *Int. J. Numer. Methods Eng.* 23 (3): 353-366. https://doi.org/10.1002/nme.1620230303.

* Shen, S.-L., H.-N. Wu, Y.-J. Cui, and Z.-Y. Yin. 2014. "Long-term settlement behaviour of metro tunnels in the soft deposits of Shanghai." *Tunnelling Underground Space Technol.* 40: 309-323. https://doi.org/10.1016/j.tust.2013.10.013.

* Shen, S.-L., Y.-X. Wu, and A. Misra. 2017. "Calculation of head difference at two sides of a cut-off barrier during excavation dewatering." *Comput. Geotech.* 91: 192-202. https://doi.org/10.1016/j.compgeo.2017.07.014.

* Sheng, D., S. Sloan, and H. Yu. 2000. "Aspects of finite element implementation of critical state models." *Comput. Mech.* 26 (2): 185-196. https://doi.org/10.1007/s004660000166.
Sołowski, W., and S. Sloan. 2015. "Evaluation of material point method for use in geotechnics." *Int. J. Numer. Anal. Methods Geomech.* 39 (7): 685-701. https://doi.org/10.1002/nag.2321.

Soundararajan, K. K. 2015. *Multi-scale multiphase modelling of granular flows*. Cambridge, UK: Univ. of Cambridge.

Staron, L., and E. Hinch. 2005. "Study of the collapse of granular columns using two-dimensional discrete-grain simulation." *J. Fluid Mech.* 545 (1): 1-27. https://doi.org/10.1017/S0022112005006415.

Utili, S., T. Zhao, and G. Houlsby. 2015. "3D DEM investigation of granular column collapse: Evaluation of debris motion and its destructive power." *Eng. Geol.* 186: 3-16. https://doi.org/10.1016/j.enggeo.2014.08.018.

Wu, Y.-X., J. S. Shen, W.-C. Chen, and T. Hino. 2017a. "Semi-analytical solution to pumping test data with barrier, wellbore storage, and partial penetration effects." *Eng. Geol.* 226: 44-51. https://doi.org/10.1016/j.enggeo.2017.05.011.

Wu, Y.-X., S.-L. Shen, and D.-J. Yuan. 2016. "Characteristics of dewatering induced drawdown curve under blocking effect of retaining wall in aquifer." *J. Hydrol.* 539: 554-566. https://doi.org/10.1016/j.jhydrol.2016.05.065.

Wu, Z.-X., Z.-Y. Yin, Y.-F. Jin, and X.-Y. Geng. 2017b. "A straightforward procedure of parameters determination for sand: A bridge from critical state based constitutive modelling to finite element analysis." *Eur. J. Environ. Civ. Eng.* https://doi.org/10.1080/19648189.2017.1353442.

Yao, Y., W. Hou, and A. Zhou. 2009. "UH model: Three-dimensional unified hardening model for overconsolidated clays." *Géotechnique* 59 (5): 451-469. https://doi.org/10.1680/geot.2007.00029.

Yao, Y., D. Lu, A. Zhou, and B. Zou. 2004. "Generalized non-linear strength theory and transformed stress space." *Sci. China Ser. E: Technol. Sci.* 47 (6): 691-709. https://doi.org/10.1360/04ye0199.

Yao, Y., D. Sun, and H. Matsuoka. 2008. "A unified constitutive model for both clay and sand with hardening parameter independent on stress path." *Comput. Geotech.* 35 (2): 210-222. https://doi.org/10.1016/j.compgeo.2007.04.003.

Yin, Z. Y., and C. S. Chang. 2013. "Stress-dilatancy behavior for sand under loading and unloading conditions." *Int. J. Numer. Anal. Methods Geomech.* 37 (8): 855-870. https://doi.org/10.1002/nag.1125.

Yin, Z.-Y., C. S. Chang, and P.-Y. Hicher. 2010. "Micromechanical modeling for effect of inherent anisotropy on cyclic behaviour of sand." *Int. J. Solids Struct.* 47 (14-15): 1933-1951. https://doi.org/10.1016/j.ijsolstr.2010.03.028.

Yin, Z.-Y., P.-Y. Hicher, C. Dano, and Y.-F. Jin. 2017a. "Modeling mechanical behavior of very coarse granular materials." *J. Eng. Mech.* 143 (1): C4016006. https://doi.org/10.1061/(ASCE)EM.1943-7889.0001059.

Yin, Z.-Y., Y.-F. Jin, S.-L. Shen, and H.-W. Huang. 2017b. "An efficient optimization method for identifying parameters of soft structured clay by an enhanced genetic algorithm and elastic-viscoplastic model." *Acta Geotech.* 12 (4): 849-867. https://doi.org/10.1007/s11440-016-0486-0.

Zenit, R. 2005. "Computer simulations of the collapse of a granular column." *Phys. Fluids* 17 (3): 031703. https://doi.org/10.1063/1.1862240.

Zhang, X., K. Krabbenhoft, D. Sheng, and W. Li. 2015. "Numerical simulation of a flow-like landslide using the particle finite element method." *Comput. Mech.* 55 (1): 167-177. https://doi.org/10.1007/s00466-014-1088-z.