Analytical and Numerical Solutions for Three-Dimensional Granular Collapses

Abstract

This research paper presents a comprehensive approach to investigating dry granular collapses in three dimensions, by combining analytical, numerical, and experimental methods. The experimental investigation utilised a novel apparatus to study granular collapses in the laboratory. It is demonstrated that a quasistatic understanding of granular collapses can accurately predict the final normalised run-out distances for dynamic granular collapses. Our approach involved establishing a correlation between the angle of repose and the initial aspect ratio of the granular column. It is also shown that the material point method (MPM) is suitable for modelling granular collapses in three dimensions. Our in-house solver was further validated using experimental evidence under an explicit formulation, resulting in good agreement between the numerical and experimental results. The findings demonstrate the effectiveness of our in-house solver for three-dimensional granular collapse modelling.

1. Introduction

The maximum angle above which noncohesive granular material starts to flow is expressed by the angle of repose ${\theta }_c$ [1,2,3]. This angle is a dominant feature of such material and is stable below the critical angle. The angle varies from 25${}^{\circ }$ for smooth particles to 45${}^{\circ }$ for angular particles [4,5,6]. Dry sands allow for typical values of approximately 35${}^{\circ }$, whereas it is much more important under wet conditions with values of approximately 90${}^{\circ }$ or even greater [7]. This angle is related to the friction coefficient $\mu$ [8], i.e., $\mu =tan\left({\theta }_c\right)$ (see [7] for further details). Other factors that influence the shape of a granular pile include gravity, grain properties (e.g., roughness, sphericity, and grain size [9,10]), the number of particles involved [11,12,13], and external perturbations such as vibrations that can cause fluidisation, leading to a relaxation stage [3].

The collapse of a dry or wet granular column [14,15,16,17] is a well-known problem, for which various experimental and numerical works [18,19,20,21] have been conducted. A fundamental metric is the initial aspect ratio of the column ${\lambda }_0=h_0/r_0$, where $h_0$ and $r_0$ are the initial height and radius of the column, respectively.

The following well-established scaling law [18,22,23,24] relates ${\lambda }_0$ to $(r_{\infty }-r_0)/r_0$, i.e.,

$$\begin{array}{c}\hfill {\displaystyle \frac{r_{\infty }-r_0}{r_0}}=\left\{\begin{array}{cc}\alpha {\lambda }_0\hfill & ,{\lambda }_0<{\lambda }_c,\hfill \\ \alpha {\lambda }_0^{\gamma }\hfill & ,{\lambda }_0\ge {\lambda }_c,\hfill \end{array}\right.\end{array}$$

(1)

where $\alpha$ and $\gamma$ are material-dependent coefficients [5,25], ${\lambda }_c$ is the initial critical aspect ratio, $r_{\infty }$ is the final radii of the column, and i.e., is the final run-out distance.

A transition phase exists, according to the initial aspect ratio, i.e., a change from a truncated cone to a cone shape [26]. References [18,19,20] proposed a power scaling law to $(r_{\infty }-r_0)/r_0$ for high ${\lambda }_0$ and linear scaling for low ${\lambda }_0$. Even though such behaviour is often observed in many studies, the value of this transition remains actively discussed.

Reference [27] proposed the following semi-empirical equation to fit their experimental results:

$$\begin{array}{c}\hfill \frac{r_{\infty }}{r_0}=\left\{\begin{array}{cc}\frac{1}{2tan\left({\theta }_c\right)}\left({\lambda }_0+{\left(4{tan}^2\left({\theta }_c\right)-\frac{{\lambda }_0^2}{3}\right)}^{1/2}\right)\hfill & ,{\lambda }_0<0.74,\hfill \\ {\left(\frac{3{\lambda }_0}{0.74}\right)}^{1/2}\hfill & ,{\lambda }_0\ge 0.74,\hfill \end{array}\right.\end{array}$$

(2)

which implies that ${\theta }_c=\mathrm{cst}.$ (constant), modulated by the variation of ${\lambda }_0$. Such a formulation was in agreement with experimental evidence, i.e., a constant final height $h_{\infty }\approx min(h_0,r_{0}tan\left({\varphi }_y\right))$, where ${\varphi }_y\approx 36.5^{\circ }$ is the internal yield angle. We refer to it as a solution to the dynamic deformations of the granular column.

This manuscript proposes an analytical solution to the final run-out distance of the granular continuum, for which we assume quasistatic deformations, i.e., $h_{\infty }=r_{\infty }tan\left({\theta }_c\right)$, whereas [27] proposed $h_{\infty }\approx min(h_0,r_{0}tan\left({\varphi }_y\right))$, i.e., dynamic deformations. These two understandings are summarised in Figure 1. In the case of quasistatic deformations, the granular pile vertically grows in proportion with ${\lambda }_0$, and we assumed that ${\theta }_c\approx {\varphi }_y$.

The analytical solution is further compared to: (i) an experimental data collection using a newly designed apparatus and (ii) the reference solution provided by Lajeunesse [27]. The expected differences between the quasistatic solution highlight the characteristic smooth transition between the quasistatic and dynamic deformations of the column. A solution to unify this quasistatic hypothesis with the dynamic understanding of [27] is then provided. Furthermore, the experimental results were compared with three-dimensional numerical solutions using the material point method (MPM), which is an effective numerical technique for modelling large deformation mechanics [28,29,30,31].

2. Materials and Methods

2.1. Laboratory Experiments

2.1.1. Experimental Setup and Data Collection

To eliminate any significant influence from the walls of the container on the flow dynamics of the granular material, a newly designed apparatus is proposed. In previous experiments, the cylinder was rapidly raised upward, resulting in uncontrolled boundary effects at the interface. The proposed experimental setup comprises a cylinder made of three independent shells that quickly open radially outward using a high-performance pneumatic system. This ensures the consistency and reproducibility of the experiments while avoiding any significant influence of the boundaries. The granular material is released and flows freely on a rough wooden surface (see Figure 2). Measurements of the final run-out distance are taken once the column has fully relaxed.

The dense front of the granular mass defines the maximum run-out distance, with measurements taken at six radial locations relative to the centre of the column. The final run-out distance is determined as the average of these six measurements, and individual grain positions are not considered as they are influenced by their previous gaseous state.

The granular mass consists of polydisperse and highly angular silicon carbide beads (SiC) with an average bead diameter of approximately $0.11$ cm and a density of $3.21$ g·cm${}^{-3}$. Assuming an equivalent spherical shape, the grain mass is estimated to be $0.02$ g. The friction coefficient is relatively high at approximately $0.77$, and it was determined by laboratory measurements outlined in the following subsection.

2.1.2. Experimental Estimation of the Angle of Repose ${\theta }_c$

The angle of repose of the granular material used in this study was measured by slowly raising a vertical cylinder initially filled with the material. The outcome of this process is the formation of a granular pile that is a valid representation of the granular material’s quasistatic relaxation. Pictures of the pile were taken using a camera, a CANON EOS 450D, and the procedure was repeated 125 times. The angle of repose over $n=250$ measurements was determined by processing the pictures and measuring the angles of both sides of the granular pile. Figure 3 shows the simple procedure to determine ${\theta }_c$.

An estimation of the maximal radial distortion was performed using the Agisoft Lens commercial software, and it indicated that the distortion at the edge of the images was approximately 0.18%, corresponding to a metric distortion of ${10}^{-4}$ mm when considering an average sensor distance of $d=600$ mm. This value was deemed sufficiently low to be negligible.

A cumulated average angle of repose was further calculated and is given by

$$\begin{array}{c}\hfill {\langle {\theta }_c\rangle }_n=\frac{1}{n}\sum _{i=1}^n{\theta }_{c,i},\end{array}$$

(3)

where $n=250$ measurements (a quantity denoted by the symbol #) of various angles of repose. An average value was identified when an equilibrium was reached, i.e., ${\partial }_n{\langle {\theta }_c\rangle }_n\to 0$. Figure 4 shows the overall results of the experimental measurements of the angle of repose. This resulted in an average angle of repose $\langle {\theta }_c\rangle =37.{55}^{\circ }\pm 0.{29}^{\circ }$. Using the relation $\mu =tan\left(\langle {\theta }_c\rangle \right)$, one obtains $\mu \approx 0.77$.

2.2. A Continuously Smooth Piecewise Analytical Solution

To define a solution $\mathcal{R}({\theta }_c,{\lambda }_0)\equiv (r_{\infty }-r_0)/r_0$ resulting from the quasistatic relaxation of a granular column as a function of $tan\left({\theta }_c\right)$ and ${\lambda }_0$, let us assume, based on Figure 5, that (i) the mass is conserved and (ii) both the angle of repose ${\theta }_c$ and the initial aspect ratio ${\lambda }_0$ govern the relaxation of the granular column, i.e., a granular pile due to quasistatic deformations during which the internal energy is dissipated by slow frictional interactions only.

The initial volume of the cylinder is given by $V_0=\pi r_0^2{h}_0$, and its final volume is given by $V_{\infty }=1/3\pi r_{\infty }^2{h}_{\infty }$ (i.e., a cone) or $V_{\infty }=1/3\pi (r_{\infty }^2+r_{\infty }\Delta r+\Delta r^2)$ (i.e., a truncated cone). As demonstrated in Figure 5, a transition occurs at ${\lambda }_c:=tan\left({\theta }_c\right)\sqrt{3}$, i.e., when ${\lambda }_0\to {\lambda }_c$.

Considering the relevant dimensions of the problem (Figure 5), the principle of mass conservation dictates that $V_i=V_{\infty }$, where it is assumed that the granular material is incompressible. The process of equating volumes and substituting and collecting terms results in the following equation:

$$\begin{array}{c}\hfill \mathcal{R}({\theta }_c,{\lambda }_0)=\left\{\begin{array}{cc}\frac{1}{2}\left(\frac{{\lambda }_0}{tan\left({\theta }_c\right)}+{\left(4-\frac{{\lambda }_0^2}{3{tan}^2\left({\theta }_c\right)}\right)}^{1/2}\right)-1\hfill & ,{\lambda }_0<tan\left({\theta }_c\right)\sqrt{3},\hfill \\ {\left(\frac{3{\lambda }_0}{tan\left({\theta }_c\right)}\right)}^{1/3}-1\hfill & ,{\lambda }_0\ge tan\left({\theta }_c\right)\sqrt{3}.\hfill \end{array}\right.\end{array}$$

(4)

This analytical solution assumes quasistatic deformations of the column. To consider the dynamic deformation of the column with a constant final height $h_{\infty }$ (see Figure 1II B), we relate the angle of repose ${\theta }_c$ to the initial aspect ratio ${\lambda }_0$, i.e., ${\theta }_c=f\left({\lambda }_0\right)$. ${\theta }_c$ is no longer constant, but depends on the initial aspect ratio of the column ${\lambda }_0$ and is given by

$${\theta }_c\left({\lambda }_0\right)={tan}^{-1}\left({\left[{tan}^3\left({\varphi }_y\right)/\left(3{\lambda }_0\right)\right]}^{1/2}\right).$$

(5)

Inserting Equation (5) within Equation (4) results in a similar formulation with respect to Equation (2) [27].

2.3. Numerical Simulation

The material point method (MPM) was originally proposed by Sulsky (1994) as an extension of the particle-in-cell method. In the MPM, the weak form of the momentum equations is solved on an Eulerian background mesh. Nodal solutions are updated and then mapped to material points, which can be considered as moving Gauss points. The state variables, such as stresses or displacements, are transported by the material points as shown in Figure 6 (taken from [32]). This makes the MPM capable of handling large deformations, such as those occurring in granular flows.

Numerical simulations were conducted to study the three-dimensional dynamics of the granular collapse. The material point method and its variant, the generalised interpolation material point method (GIMP), were employed [33]. An explicit MPM solver that takes advantage of modern graphics processing unit (GPU) architectures was implemented, i.e., ep2-3De v1.0 (The latest version of the solver is available for download from GitHub at: https://GitHub.com/ewyser/ep2-3De (accessed on 10 August 2021)). Further details of the implementation of the solver can be found in [32,34]. Because of the large deformations involved during the collapse, we selected the uGIMP variant, i.e., the material point’s spatial extent is constant (i.e., undeformed; see [35] for the limitations of this variant). A non-associated Drucker–Prager plastic rheology was selected to reproduce the granular collapse [36].

2.4. Numerical Parameters and Geometry

Recently, Reference [37] demonstrated that density and stiffness properties have a negligible effect on the morphology and run-out distance of granular collapses. Therefore, the same values as [37] were used as we considered the material properties for an effective continuum medium, i.e., a density $\rho =2000$ kg·m${}^{-3}$, with Young’s modulus $E=5.84$ MPa, a cohesion $c=0$ kPa, and a Poisson’s ratio $\nu =0.3$. The friction angle was $\varphi =37.{55}^{\circ }$, i.e., $\varphi =\langle {\theta }_c\rangle$, which was measured from our experiments. A local damping D was introduced, for which $D\in [0.05;0.1]$ is a commonly accepted range for explicit formulations [38]. This local damping is proportional to the magnitude of the out-of-balance forces calculated on the background mesh [38].

Even though GPU programming enables performant (in terms of wall-clock time) numerical solvers and one of the major limitations is the required memory [34], which can far exceed the hardware limit, only one quarter of the granular column was considered, assuming that the horizontal momentum transfers are sufficiently small to be neglected. This allowed an important amount of memory utilisation to be spared on the GPU during computation.

Since the initial aspect ratio ${\lambda }_0$ of the column strongly governs the run-out, we assumed that the numerical geometry could differ from the experimental setting. Consequently, a column of radius $r_0=1$ m was considered, to artificially increase the numerical time steps, which are restricted by the Courant–Friedrichs–Lewy (CFL) condition, i.e., adaptive time steps are implemented in ep2-3De v1.0. The three-dimensional background mesh is made of regular quadrilateral elements. Roller boundary conditions, i.e., free-slip boundary conditions, were enforced on the side boundaries of the background mesh, whereas a no-slip boundary condition was enforced on the bottom. The granular column was discretised by 40 elements along the x and y directions, and $n_{pe}=8$ were regularly assigned per the initially filled element (see

Table 1. The granular column was discretised by 40 elements along the x and y directions, and $n_{pe}=8$ were assigned per the initially filled element.

${\mathit{\lambda }}_0$	${\mathit{n}}_{\mathit{el}}$	${\mathit{n}}_{\mathit{no}}$	${\mathit{n}}_{\mathit{mp}}$	$\mathit{\Delta }\mathit{x}$ [m]
0.5	27,104	30,375	25,120	0.05
1.0	169,344	180,625	50,240	0.05
2.0	1,183,424	1,225,125	100,480	0.05
4.0	8,817,984	8,978,125	200,960	0.05

). The background mesh depends on the initial height of the column, and it is defined to be sufficiently large to fully enclose the collapse without the influences of side boundaries.

3. Results

3.1. Analytical Solutions and Experimental Collapses

Figure 7a shows the experimental results from the laboratory experiments fit by: (i) the analytical solution under the quasistatic deformation hypothesis of the column and (ii) the solution given by [27]. The latter agrees with the experimental granular collapses, whereas the former (the present quasistatic understanding of the deformation) rapidly diverges from the experimental data.

In this case, ${\theta }_c={\varphi }_y$ was considered, where the internal yield angle was ${\varphi }_y={tan}^{-1}\left(\mu \right)$ with $\mu =0.77$. The latter was inferred with the experimental protocol previously presented. When considering ${\theta }_c=f\left({\lambda }_0\right)$, the proposed analytical solution (green line in Figure 7a) is then in agreement with both the experimental results and Equation (2) [27].

Regarding the experimental data, a small amplitude of the error measurements was observed (see Figure 7b). This indicates that the present observations are consistent and reliable.

To quantify the goodness of fit of Equation (4), a power-law fit was applied to the experimental data. It is given by

$$\begin{array}{c}\hfill {\displaystyle \frac{(r_{\infty }-r_0)}{r_0}}=0.63{\lambda }_0^{1.30},\end{array}$$

(6)

for which the root-mean-squared error $\mathrm{RMSE}=0.13$ ($R^2=0.99$), which is lower than the proposed analytical solution. An $\mathrm{RMSE}=0.23$ was obtained for ${\varphi }_y={tan}^{-1}\left(\mu \right)$. The optimal ${\varphi }_y$ was further investigated (see Figure 8) for which the $\mathrm{RMSE}$ for Equation (4) was the smallest when fit to the experimental data.

The difference observed for the parameter $\mu$ (i.e., 0.71 and 0.77 according to Figure 8) appears reasonable. The value inferred from the robust fitting showed the smallest RMSE; however, the value inferred from the laboratory measurements has a physical meaning (i.e., the friction) despite its greater error. Since the error difference is acceptable, we considered $\mu =0.77$.

3.2. Experimental and Numerical Granular Collapses

A typical numerical solution obtained with the uGIMP variant in ep2-3De v1.0 (see [34]) is shown in Figure 9. The initial aspect ratio was ${\lambda }_0=2$, and the local damping was $D=0.05$. The granular mass spread in a realistic fashion on the bottom surface. The equivalent plastic strain ${ϵ}_{\mathrm{eqv}}^p$ highlighted intense zones of shearing. However, successive shear bands were roughly resolved because of the numerical resolution, i.e., $\Delta x=0.05$ m. One can observe that most of the plastic deformation was superficial. In addition to the overall elastoplastic deformations of the granular column, shallower granular avalanches also occurred. Few material points can be considered fully disconnected from the main body, i.e., the gaseous state observed during the experiments. Similarly, the proportion of material points in a disconnected state increased as the local damping was reduced.

Different values were selected for the local damping, i.e., $D=\{0.0,0.05,0.1\}$ (see Appendix A for further details). For simplicity, the maximal radial distance of the farthest material point is reported to determine $r_{\infty }$. For $D=0.05$ (see Figure 10), the numerical solution agreed well with the experimental granular collapses. If the damping was smaller or greater than 0.05, the numerical model either overestimated or underestimated the normalised run-out distance, respectively. As such, local damping is an important parameter and requires an iterative calibration process.

Increasing the local damping even more (i.e., $D=0.4$) led to a numerical solution closer to the proposed analytical solution assuming quasistatic deformation. This makes sense since most of the out-of-balance forces were damped out during the calculation, yielding final run-out distances close to the quasistatic state assumed previously.

4. Discussion

The difference between the present analytical solutions (quasistatic and dynamic hypotheses; see Figure 11) expresses an important principle under the following hypothesis. The kinetic energy loss nearly asymptotically increased during the collapse because of an increase in the elastoplastic collision rate. This suggested that the final angle of repose ${\theta }_c$ of a granular collapse expresses the amount of energy lost during the process, compared to a purely quasistatic deformation of the column: the greater the difference between ${\theta }_c$ and ${\varphi }_y$, the greater the kinetic energy loss is.

Further interpretation yielded the following. The near-asymptotic behaviour of $tan\left({\theta }_c\right)$ (see Figure 11) could be understood because of a microscopic steady-state collision rate during the granular collapse. The collision rate rapidly increased as the initial aspect ratio increased. When the grains reached their free-fall velocity, their collision rate became steady. When such an equilibrium is resolved, the energy dissipation rate within the system can no longer change.

Within the MPM framework under an explicit formulation, the numerical solutions agreed with the experimental results of granular collapses. This confirmed the solver ep2-3De v1.0 to be an appropriate tool for the numerical modelling of dynamic granular collapses. However, a proper calibration procedure for local damping must be conducted. An increase of the local damping yielded a numerical solution closer to the quasistatic analytical solution. However, the explicit formulation of the solver is not well suited to further investigate such numerical transitions. An implicit formulation [39,40] should be preferred to fully resolve quasistatic granular collapses.

Because of hardware limitations (e.g., on-chip memory limit), a full three-dimensional model of granular collapse is not yet possible. This should be the focus of future studies by proposing a multi-GPU implementation of the solver ep2-3De v1.0 using a message-passing interface standard (MPI), such as Open MPI.

As suggested by [37], some material properties (i.e., stiffness and density) had little influence on the behaviour of the collapse. The only common material parameter between the experiments and the numerical model was the friction angle. Nevertheless, the numerical solutions agreed well with the experimental data. This also demonstrated that the only geometrical parameter that truly matters was the initial aspect ratio ${\lambda }_0$ of the column.

5. Conclusions

We proposed an analytical solution for the normalised run-out distance of three-dimensional quasistatic granular collapses. We further introduced a correction to consider dynamic collapses and the influence of the initial aspect ratio ${\lambda }_0$ of the column as a modulation of the angle of repose ${\theta }_c$. Such an analytical solution was found to be in good agreement with the experimental results. This demonstrated that a quasistatic understanding of granular collapses can include its dynamic counterpart as well by including a relation between the final angle of repose with the initial aspect ratio of the column. Here, we can only further assume that energy dissipation during the collapse plays a key role in determining the final angle of repose. This may explain the observed variations in the final angle of repose.

Furthermore, we validated our proprietary explicit solver, ep2-3De v1.0, by comparing it with experimental granular collapses. Our results demonstrated good agreement and revealed that the introduction of a damping factor modulation enabled us to simulate both quasistatic and dynamic granular collapses. This also allowed us to establish a calibration procedure to determine the optimal damping parameter for use within an explicit material point framework. Taken together, these findings provide strong evidence for the accuracy of our solver in resolving granular collapses. However, future work should focus on two important directions: (1) implementing a multi-GPU approach to overcome the limitations of on-chip memory and (2) developing an implicit formulation of the solver that can natively resolve quasistatic granular collapses and facilitate comparisons with damped explicit solutions.

To improve the robustness and applicability of the proposed analytical solution, future experimental investigations should explore the effects of various shapes of grains. Investigating the impact of irregular grain shapes on the validity of the proposed analytical solutions would be particularly interesting. This objective could be accomplished through the use of different types of granular materials, e.g., long and elongated grains such as rice or quite spherical beads.