Modeling sand slides by a mechanics-based degenerate parabolic equation

Abstract

Avalanching plays a crucial role in granular materials dynamics, in particular in the evolution of the shape of the leeward side of sand dunes. This paper presents a physically-based mathematical model capable of reproducing the kinematic evolution of the surface of sand piles and to obtain eventually the stationary configurations, in the presence of external sources as well. Simulation results with different boundary conditions and geometries are reported in order to show the high flexibility of the model. The model is also validated by means of comparison with the experimental results of different authors.

1. Introduction

One of the crucial phenomena that contributes to the evolution of sand piles is the formation of avalanches in regions with steep slopes. In fact, for instance, the sand that is transported by the wind in dunes accumulates and then spontaneously slips when the steepest slope exceeds the angle of repose, which in general depends on the properties of the media, firstly grain size and density, and in the case of dry sand it is close to 32°. Experimental approaches are useful to obtain an estimation of the angle of repose, but they are hard to be used in complex processes involving the interaction of wind with the substratum, for instance, to evaluate the evolution of the sand bed level, or of the dune eventually interacting with other structures.

Such mathematical models are also of interest for other similar natural phenomena, such as snow avalanches, landslides, debris accumulation, scour and sediment transport, and also in industrial applications, for instance in the stocking processes of many different materials including rocks, seeds, and all media composed of discrete particles.

Our interest in understanding and simulating the dynamics of sand slides is part of a wider effort consisting in deducing a complete model that also includes sand erosion, transport, and sedimentation [1]. In such a model, sand deposition and avalanches together with sand erosion contribute to determining the free boundary between the sand bed and the multiphase flow composed of sand grains transported by the wind, which includes the saltation layer in the region closer to the sand bed. In this framework, the specific need addressed in this article is to deduce a sand avalanche model that is able to describe the evolution of the sand bed in a simultaneously accurate, simple, and fast way. The final engineering aim of the complete model is then to devise and numerically test the efficacy of protection measures for human infrastructure in arid environments, such as roads and railways, pipelines, industrial facilities, farmland, towns, and buildings.

Detailed avalanche models have been developed in the past, starting from the pioneering papers by Savage and Hutter [2], in which hydrodynamical models, based on Saint-Venant equations, were originally proposed. They started from the incompressibility condition and momentum conservation equations in the flowing layer and then integrated them over the thickness of the rolling layer, not considering erosion and deposition. This gap was then filled, for instance, by Douady et al. [3], Khakhar et al. [4], and Gray [5]. A detailed analytical study of Savage and Hutter’s model has been performed by Colombo et al. [6]

A variational approach has also been suggested. Aronsson et al. [7] introduced the p -Laplacian operator in order to model a non-homogeneous superficial diffusion. Prigozhin [8] developed a variational model based on mass balance for the sand pile and assuming the surface flow is directed down the path of steepest descent subject to a slope constraint. Prigozhin and Zaltzman [9] allowed the coefficient in front of the gradient of the sand pile to vanish if the slope is below a threshold. Other successful approaches are based on the theory of self-organized criticality [10] and cellular automata [11].

Bouchaud et al. [12] paved the way to the so-called two-layer models. This modeling approach starts from the assumption that sand piles can be divided into a static and a rolling layer and include source terms as well. The model consequently uses two variables, the local height of the static layer and the local density (or height) of the rolling layer, linked by an exchange term that allows static grains to start moving. The drift downhill velocity of the rolling grains is considered constant. A modified version of the two-layer model has been proposed by de Gennes and coauthors [13–16]. Prigozhin and Zaltzman [9] overcame the constant velocity hypothesis by assuming a proportionality with the gradient of the height of the sand pile. In the models proposed by Hadeler and Kuttler [17, 18], velocity is usually constant, though it is mentioned that it might depend on the slope (see, for instance, [19, 20]). An exchange term proportional to the difference between the slope of the sand bed and the angle of repose governs the transfer of sand between the static and the rolling layer.

A detailed analytical study of this type of models has been performed by Cannarsa and coworkers [21–24] for different boundary conditions, whereas Falcone and Finzi Vita [25–27] also focused on the proper numerical integration method.

In spite of the fact that two-layer models can describe several situations well, they present some drawbacks such as the presence of artifacts when the boundary of the domain changes type (from a wall to an open end) [23, 26], the difficulty in controlling what happens at the boundary due to the hyperbolic character of the equations, and the formation of artificial valleys on the top of the sand pile due to the discontinuity of downhill velocities on the two slopes, or in situations in which the maximum height of the sand pile should be reached close to a wall [28].

A further important problem occurs when applying the model to configurations with supercritical initial profiles, that is when the sand pile initially has a portion of the slope with an angle larger than the angle of repose, as in the experiments in [29–36], where the lateral surfaces of the sand pile are vertical and there is no external sand flux. In this case, in general, it is hard to identify the proper initial conditions to be joined to the two-layer models. Moreover, it can be shown that the evolution and the final configuration crucially depend on how the initial profile is divided into a static and a rolling part. Paradoxically, if it is assumed that at the beginning all the sand is static, as it would look natural to do, then the system would not evolve at all.

In this work, we deduce a simple model, that for sake of convenience we will call Degenerate Parabolic Sliding Model and shorten as DPSM, for the height of the sand pile deduced on the basis of a mass balance equation under the assumption that the thickness of the creep layer is small and that sand grains in it flow down the direction of the steepest descent due to gravity. We show that the stationary configurations satisfy the eikonal equation, as classically known and previously found in other two-layer models (see [17]). The model is flexible with respect to boundary conditions and influx terms and is experimentally validated against column collapse results in [29, 30, 35, 36].

This paper is organized as follows. In Section 2, we derive the model based on a suitable physically-based closure of the mass balance equation. Discretization and numerical aspects are dealt with in Section 3, where the meaning of the model is also described from a discrete point of view. Then the results from the numerical simulation of the degenerate parabolic sliding model are reported and validated against several experimental results in Section 4. In Section 5 we draw some conclusions.

2. The degenerate parabolic sliding model

Consider a control volume $V$ with a vertical lateral surface $\partial V_{\ell }$ , the basis $A$ deep in the sand (where the sand does not move) and the top $A_t$ high in the air. Referring to Figure 1, the volume is divided by the surface $z=h\left(x,y,t\right)$ into two regions, the sand bed $V_s$ and the air above it $V_a$ . The sand bed is further divided into a lower zone in which the sand does not move and has a close packing density ${\rho }_{*}$ and a thin creep layer $S$ of thickness $\delta$ where sand can move and has a density $\rho \left(x,y,z,t\right)$ that is close to the close packing one. Hence, the mass of sand in $V_s$ is

$$M={\int }_V\rho \mathrm{dV}={\int }_A\left\{{\rho }_{*}\left[h-\frac{\delta }{|\mathbf{k}\cdot \mathbf{n}|}\right]+{\int }_{h-\frac{\delta }{|\mathbf{k}\cdot \mathbf{n}|}}^h\rho \mathrm{dz}\right\}\mathrm{dA},$$

where $\mathbf{n}$ is the normal to the surface and $\mathbf{k}$ is the vertical unit vector. Using the mean value theorem for integrals we have that

$$M={\int }_A\left\{{\rho }_{*}h\left(\mathbf{x},t\right)-\left[{\rho }_{*}-\rho \left(x,y,ẑ,t\right)\right]\frac{\delta }{|\mathbf{k}\cdot \mathbf{n}\left(\mathbf{x},t\right)|}\right\}\mathrm{dA},$$

where $ẑ\in \left[h\left(\mathbf{x},t\right)-\frac{\delta }{|\mathbf{k}\cdot \mathbf{n}\left(\mathbf{x},t\right)|},h\left(\mathbf{x},t\right)\right]$ . Because of both the smallness of $\delta$ and of the closeness of the density in the creep layer to the close packing one, the second term in parenthesis can be simplified and one has that

$$M={\int }_{\partial V_{\ell }\cap \partial V_a}{\rho }_{*}h\left(\mathbf{x},t\right)\mathrm{dA}.$$

Figure 1.

Definition of the domains and the subdomains $V_a$ colored in pale blue and $V_s$ colored in darker and ligther orange (corresponding to the creep layer).

Figure 2.

Evolution of the sand pile profile on an open table due to a constant source q = 0.05 located in the sub-intervals $D_w$ highlighted by the dashed segments in the upper part of the graphs. Specifically: (a) $D_w=$ (0, 1); (b) $D_w=$ (0, 0.4); (c) $D_w=$ (0.2, 0.4) ∪ (0.8,1); (d) $D_w=$ (0.18, 0.22) ∪ (0.78,0.82). Stationary profile is shown in red.

Even if in the following we will assume that some sand may be transported through the boundary of $V_a$ we will also assume that the mass that enters through it readily sediments, so that the air can always be assumed to be clear. This assumption has no effect on the final model.

Balance of mass states that the mass of the sand contained in the volume can change due to the flux of sand through the boundary. In order to evaluate it we distinguish between the flux in the creep layer and that coming through the air, which can be written as

$$\begin{array}{c}{\int }_{\mathrm{\partial V}\cap \partial V_a}\mathbf{q}(\mathbf{x},t)\cdot \mathbf{n}d\Sigma ={\int }_{A_t}\mathbf{q}(x,y,z=H,t)\cdot \mathbf{k}\mathrm{dA}\\ +{\int }_{\partial V\cap \partial V_a}\mathbf{q}(\mathbf{x},t)\cdot \mathbf{n}d\Sigma ,\end{array}$$

where the first term is the outflux through the top of the cylinder and the second is that through the lateral surface.

The main hypothesis to evaluate the contribution due to the motion of the sand grains in the creep layer S is that sand grains move in the direction $\mathbf{t}$ of the steepest descent on the surface (see the Appendix for its computation), i.e.,

$${\int }_{\partial V_{\ell }\cap \partial V_s}\mathbf{q}\cdot \mathbf{n}d\Sigma ={\int }_{\partial V\cap \partial S}\rho w\mathbf{t}\cdot \mathbf{n}d\Sigma ,$$

(1)

where the speed w will be then identified as the limit velocity of a body sliding down a slope under the action of gravity, friction, and a drag force taken to be proportional to its speed, due to the impact with the other grains of the sand bed. This means that

$$\mathbf{v}=w\mathbf{t},\mathrm{with}\mathbf{t}=-\frac{\nabla h+|\nabla h{|}^2\mathbf{k}}{|\nabla h|\sqrt{1+|\nabla h{|}^2}},$$

(2)

Coming to the sliding speed w, simplifying consistently the microscopic structure of the sand bed, by a simple body diagram on the single sand grain we assume that, on average, motion is triggered if the component of the gravitational force along the surface ( $\mathrm{mg}\sin \theta$ where θ is the local angle of the steepest slope of the surface and m is the mass of the grain) is larger than the maximum static friction force ( ${\mu }_s\mathrm{mg}\cos \theta$ where ${\mu }_s$ is the friction coefficient). Under such conditions and hypotheses one has that the limit sliding speed is identified by $-\gamma w+\mathrm{mg}\sin \theta -{\mu }_d\mathrm{mg}\cos \theta =0$ where ${\mu }_d$ is the dynamic friction coefficient. However, looking at the scale of the sand grains, the friction angles are not so well defined, which leads to some uncertainty in their experimental evaluation and to a smoother transition from rest to the limit speed for slopes closer to the static angle. For sake of simplicity, in the following we will assume that the two friction coefficients are equal and then define the limit speed as

$$w=W\cos \theta {(\tan \theta -\tan \bar{\theta })}_{+},$$

(3)

where $\bar{\theta }$ is the angle of repose, and ${(\cdot )}_{+}=\frac{(\cdot )+|\cdot |}{2}$ is the positive part of the argument in parentheses.

Equation (3) states, for instance, that if the local slope is subcritical (i.e. with an angle that is smaller than the angle of repose), then sand grains will not move. In the supercritical case (i.e. when the angle of the slope is larger than the angle of repose) the speed increases with the slope. The speed increases when the energy loss mechanisms due to the interactions among the grains, modeled through w, decreases.

An alternative way to obtain the sliding speed is to describe the motion in the creep layer as the flow of a Bingham fluid. In this case the bilinear behaviour expressed in Eq. (3) is replaced by a parabolic behaviour above a threshold value related to the yield stress. More in general, using a Herschel-Bulkley fluid with a power γ would lead to a relation with a positive part term elevated to the power $1+\frac{1}{\gamma }$ . Such closure possibilities will be studied and compared in a future work.

Taking into account of the relationship between θ and $\nabla h$ , Equation (3) can be rewritten as

$$w=W\frac{{\left(\right|\nabla h|-\tan \bar{\theta })}_{+}}{\sqrt{1+|\nabla h{|}^2}}.$$

(4)

Computing then the flux contribution (1), one has

$$\begin{array}{ccc}{\int }_{\partial V\cap S}\rho w\mathbf{t}\cdot \mathbf{n}d\Sigma & =& {\int }_{\partial A}{\rho }_{*}w\mathbf{t}\cdot \mathbf{n}\frac{\delta }{\cos \theta }d\mathrm{\Gamma ,}\end{array}$$

(5)

where using the same argument used for the evaluation of the sand flux, the density of the grains was approximated with the close packing density ${\rho }_{*}$ . In writing Eq. (5) we also used the fact that the thickness of the creep layer (normal to the sliding surface) is assumed constant and then the stripe on the lateral surface of the cylindrical control volume has height δ∕cosθ. However, taking into account that the normal to the lateral surface of the control volume is in the horizontal plane, one has that

$$\mathbf{t}\cdot \mathbf{n}=-\frac{1}{\sqrt{1+|\nabla h{|}^2}}\frac{\nabla h}{|\nabla h|}\cdot \mathbf{n},$$

and, therefore,

$$\begin{array}{c}{\int }_{\partial V\cap \partial S}\mathbf{v}\cdot \mathbf{n}d\Sigma =-{\int }_{\partial A}w\frac{\delta }{\cos \theta }\frac{1}{\sqrt{1+|\nabla h{|}^2}}\frac{\nabla h}{|\nabla h|}\cdot \mathbf{n}d\Gamma \\ =-{\int }_A\nabla \cdot \left(w\frac{\delta }{\cos \theta }\frac{1}{\sqrt{1+|\nabla h{|}^2}}\frac{\nabla h}{|\nabla h|}\right)\mathrm{dA}.\end{array}$$

Mass balance then leads to

$${\int }_A\frac{\partial h}{\partial t}\mathrm{dA}=-{\int }_{A_t}\frac{\mathbf{q}\cdot \mathbf{k}}{{\rho }_{*}}\mathrm{dA}-{\int }_{\partial V_{\ell }\cap \partial V_a}\frac{\mathbf{q}\cdot \mathbf{n}}{{\rho }_{*}}d\Sigma +{\int }_A\nabla \cdot \left(w\delta \frac{\nabla h}{|\nabla h|}\right)\mathrm{dA},$$

that by the usual continuum mechanics methods can be localized to

$$\frac{\partial h}{\partial t}=\nu \nabla \cdot \left[\frac{{\left(\right|\nabla h|-\tan \bar{\theta })}_{+}}{\sqrt{1+|\nabla h{|}^2}}\frac{\nabla h}{|\nabla h|}\right]+q,$$

(6)

where ν = Wδ and $q=-\mathbf{q}\cdot \mathbf{k}∕{\rho }_{*}$ is the net flux between possible external mass supply due to sand sedimentation in the air and mass leaving the control volume due to erosion.

The evolution equation for the height of the sand pile is then a parabolic equation that may degenerate if the maximum slope in a point forms an angle with the horizontal plane that is smaller than the angle of repose. As done before, in the following we will refer to this point (or subregion) as subcritical, whereas the complementary situation will be called supercritical. In the following, we will also extend these definitions to the entire sand pile, calling subcritical distribution of sand a configuration in which every point is subcritical. In this case, the evolution equation reduces simply to

$$\frac{\partial h}{\partial t}=q,\mathrm{for}(\mathbf{x},t):|\nabla h(\mathbf{x},t\left)\right|\le \tan \bar{\theta }.$$

(7)

We use the term supercritical distributions of sand for those cases characterized by the presence of regions where the situation is supercritical. In supercritical regions, one has the parabolic equation

$$\frac{\partial h}{\partial t}=\nu \nabla \cdot \left[\frac{|\nabla h|-\tan \bar{\theta }}{\sqrt{1+|\nabla h{|}^2}}\frac{\nabla h}{|\nabla h|}\right]+q,\mathrm{for}(\mathbf{x},t):|\nabla h(\mathbf{x},t\left)\right|>\tan \bar{\theta },$$

(8)

degenerating on the border of the time-dependent subdomain.

Looking for stationary configurations, it can be trivially observed that any $h\left(\mathbf{x}\right)$ with $|\nabla h|\le \tan \bar{\theta }$ satisfies (6) and (7) for q=0 and is then a candidate to be a stationary solution. Hence, there is non-uniqueness of the static solutions that is then identified by the evolution of the system. This is not a surprise from the physical point of view, because as is well-known, sand can be distributed in an infinite number of ways as long as the slope does not exceed the critical one.

One can then define the critical (or limit) stationary solution as the one satisfying exactly the eikonal equation

$$|\nabla h|=\tan \bar{\theta },$$

everywhere. It is well known (see, for instance, Hadeler and Kuttler [17]) that the solutions of the eikonal equations are related to the shape of the larger sand pile that can be obtained by pouring sand on an open plate. Thus, for instance, one has a cone if the plate is a circle and a pyramid if the plate is a regular polygon.

3. Finite difference discretization

In this section, we discretize Equation (6) and we will discuss the resulting discrete model showing that it has a specific physical interpretation. For the sake of clarity, we consider the one-dimensional version of Equation (6) without an external source

$$\frac{\partial h}{\partial t}=\nu \frac{\partial }{\partial x}\left[\frac{{\left(\left|\frac{\partial h}{\partial x}\right|-\tan \bar{\theta }\right)}_{+}}{\sqrt{1+{\left|\partial h\partial x\right|}^2}}\mathrm{sgn}\left(\frac{\partial h}{\partial x}\right)\right],$$

(9)

where $\mathrm{sgn}(\cdot )$ is the signum function, as the general case is a straightforward generalization.

The last term inside the square brackets gives the direction of flux and can be replaced by

$$\frac{\partial h}{\partial t}=\nu \frac{\partial }{\partial x}\left[\frac{{\left(\frac{\partial h}{\partial x}-\tan \bar{\theta }\right)}_{+}}{\sqrt{1+{\left|\frac{\partial h}{\partial x}\right|}^2}}-\frac{{\left(-\frac{\partial h}{\partial x}-\tan \bar{\theta }\right)}_{+}}{\sqrt{1+{\left|\frac{\partial h}{\partial x}\right|}^2}}\right].$$

(10)

Considering a uniform spatial discretization, calling $\alpha =\tan \bar{\theta }$ , the right-hand side of Equation (10) reads

$$\begin{array}{c}\frac{\nu }{\Delta x}\left[\frac{{\left({\left|\frac{\partial h}{\partial x}\right|}_{x_{i+1}}-\alpha \right)}_{+}}{\sqrt{1+{\left|\frac{\partial h}{\partial x}\right|}_{x_{i+1}}^2}}\frac{\frac{\partial h}{\partial x}{|}_{x_{i+1}}}{{\left|\frac{\partial h}{\partial x}\right|}_{x_{i+1}}}-\frac{{\left({\left|\frac{\partial h}{\partial x}\right|}_i-\alpha \right)}_{+}}{\sqrt{1+{\left|\frac{\partial h}{\partial x}\right|}_{x_i}^2}}\frac{\frac{\partial h}{\partial x}{|}_{x_i}}{{\left|\frac{\partial h}{\partial x}\right|}_{x_i}}\right]\\ =\frac{\nu }{\Delta x}\left[\frac{{\left(\frac{h_{i+1}^n-h_i^n}{\Delta x}-\alpha \right)}_{+}}{\sqrt{1+{\left|\frac{h_{i+1}^n-h_i^n}{\Delta x}\right|}^2}}-\frac{{\left(-\frac{h_{i+1}^n-h_i^n}{\Delta x}-\alpha \right)}_{+}}{\sqrt{1+{\left|\frac{h_{i+1}^n-h_i^n}{\Delta x}\right|}^2}}\right.\\ \left.-\left(\frac{{\left(\frac{h_i^n-h_{i-1}^n}{\Delta x}-\alpha \right)}_{+}}{\sqrt{1+{\left|\frac{h_i^n-h_{i-1}^n}{\Delta x}\right|}^2}}-\frac{{\left(-\frac{h_i^n-h_{i-1}^n}{\Delta x}-\alpha \right)}_{+}}{\sqrt{1+{\left|\frac{h_i^n-h_{i-1}^n}{\Delta x}\right|}^2}}\right)\right].\end{array}$$

Rearranging terms and introducing time discretization, Equation (10) becomes

$$\begin{array}{c}\frac{h_i^{n+1}-h_i^n}{\Delta t}=\frac{\nu }{\Delta x}\left[\frac{{(\Delta h_{i+1}^n-\Delta x\alpha )}_{+}}{\sqrt{\Delta x^2+{(\Delta h_{i+1}^n)}^2}}\right.\\ -\frac{{(-\Delta h_{i+1}^n-\Delta x\alpha )}_{+}}{\sqrt{\Delta x^2+{(\Delta h_{i+1}^n)}^2}}\\ +\frac{{(-\Delta h_i^n-\Delta x\alpha )}_{+}}{\sqrt{\Delta x^2+{(\Delta h_i^n)}^2}}\\ \left.-\frac{{(\Delta h_i^n-\Delta x\alpha )}_{+}}{\sqrt{\Delta x^2+{(\Delta h_i^n)}^2}}\right]\end{array}$$

(11)

where $\Delta h_{\alpha }^n=h_{\alpha }^n-h_{\alpha -1}^n$ .

In the general case, in which the i -th node has $N_n^i$ neighbors, Equation (11) reads

$$\frac{h_i^{n+1}-h_i^n}{\Delta t}=\nu \sum _{j=1}^{N_n^i}\frac{1}{d_{\mathrm{ij}}}\frac{{\left(\right|h_j^n-h_i^n|-d_{\mathrm{ij}}\alpha )}_{+}}{\sqrt{d_{\mathrm{ij}}^2+{(h_j^n-h_i^n)}^2}}\mathrm{sgn}(h_j^n-h_i^n),$$

(12)

where $d_{\mathrm{ij}}$ is the distance between the two nodes i and j. The stationary solution of this discrete equations can be seen as the solution of a discrete model for avalanching, constructed requiring only mass conservation and angle of repose condition compliance. Assuming unit density and unitary base area, the mass associated with a node is simply equal to the value of the height $h_i$ . Therefore, considering two generic nodes i and j, mass conservation and the angle of repose condition read

$$\left\{\begin{array}{c}{h}_i^{n+1}+h_j^{n+1}=h_i^n+h_j^n,\hfill \\ h_i^{n+1}-h_j^{n+1}=\pm d_{\mathrm{ij}}\tan \bar{\theta },\hfill \end{array}\right.$$

(13)

whose solution is

$$\left\{\begin{array}{c}{h}_i^{n+1}=\frac{h_i^n+h_j^n\pm d_{\mathrm{ij}}\tan \bar{\theta }}{2},\hfill \\ h_j^{n+1}=\frac{h_i^n+h_j^n\mp d_{\mathrm{ij}}\tan \bar{\theta }}{2},\hfill \end{array}\right.$$

(14)

which can be seen as the stationary solution of Equation (12) for the couple i –j.

In next section we apply this discretization scheme using a uniform Cartesian grid. Each node can exchange mass with all the nodes sharing with it a face of the grid. Actually, the scheme can be used with others grid as well, as it is based on mass exchange between neighbours. For stabilization, the mass exchange allowed in each iterative step is bounded as described in [37].

4. Numerical simulations

In this section, the degenerate parabolic model (6) is applied in one- and two-dimensional domains with different boundary conditions in order to show its robustness. Two types of boundaries can be conceived: an open end $\partial D_{\mathrm{open}}$ where sand can freely run away and a closed end $\partial D_{\mathrm{wall}}$ where sand needs to stop. The former refers, for instance, to the case in which sand is poured onto a table and as soon as it reaches the end of the table, it drops down. In this case the proper boundary condition to impose on Equation (6) is

$$h(\mathbf{x},t)=0,\mathrm{for}\mathbf{x}\in \partial D_{\mathrm{open}}.$$

The latter refers to the case in which the domain ends with a wall. In this case, the proper boundary condition is

$$\mathbf{v}(\mathbf{x},t)\cdot \mathbf{n}=0,\mathrm{for}\mathbf{x}\in \partial D_{\mathrm{wall}},$$

or

$${\left(\right|\nabla h(\mathbf{x},t)|-\tan \bar{\theta })}_{+}\nabla h(\mathbf{x},t)\cdot \mathbf{n}=0,\mathrm{for}\mathbf{x}\in \partial D_{\mathrm{wall}}.$$

This boundary condition can be satisfied either when $|\nabla h|\le \tan \bar{\theta }$ , or when ∇h ⋅n = 0, but in the first case Equation (6) reduces to Equation (7) that does not really call for a boundary condition. In conclusion, in the case of a wall the proper boundary condition to join to (6) is

$$\mathbf{n}\cdot \nabla h(\mathbf{x},t)=0,\mathrm{if}|\nabla h(\mathbf{x},t\left)\right|>\tan \bar{\theta },\mathrm{for}\mathbf{x}\in \partial D_{\mathrm{wall}}.$$

In the following, two main problems will be examined in detail. In the first, discussed in this section, mass is supplied in some subregions of a domain that is surrounded either by open boundaries, or by walls, or by a combination of them. In the second problem, discussed in the following section, evolution starts from a pile of sand that is supercritical somewhere and there is no mass supply. Specifically, we will start from a cylindrical pile of sand redistributing its mass on an infinite plane in order to compare the outcome of the degenerate parabolic model with existing experimental results.

Before going to the simulations, in order to understand the structure of the expected solution, let us consider the case in which the mass supply q is constant everywhere. If there are walls all around the domain, then the solution is trivially a constant h rising at a velocity q. In the opposite case, if all the borders are open, one can look for an approximate stationary solution when q ≪ ν∕L where L is the characteristic dimension of the domain. Consider first the one-dimensional problem in the interval [0,L]. The solution vanishes on the border x=0 and then increases. This suggests we should look for a solution in [0,L] of the type $h(x,t)=\tan \bar{\theta }x+ϵ\eta (x,t)$ where ϵ = qL∕ν. In this case the partial differential equation (PDE) linearizes to

$$\frac{\partial }{\partial \tilde{x}}\left(\cos \bar{\theta }\frac{\partial \tilde{\eta }}{\partial \tilde{x}}\right)+1=0,$$

where distances are made dimensionless with L. Hence, in dimensional terms,

$$h(x,t)=\tan \bar{\theta }x+\frac{q}{2\nu \cos \bar{\theta }}x(L-x),\mathrm{for}x\in \left[0,\frac{L}{2}\right],$$

(15)

and symmetrically

$$h(x,t)=\tan \bar{\theta }(L-x)+\frac{q}{2\nu \cos \bar{\theta }}x(L-x),\mathrm{for}x\in \left[\frac{L}{2},L\right].$$

(16)

Hence, the steady state is locally concave with an almost triangular shape.

The above discussion also gives an idea of how the solution tends to this stationary state (see Figure 3(a)). In fact, in the central region the sand pile grows with a flat profile rising with velocity q, then growing in an almost trapezoidal shape

$$h(x,t)=\left\{\begin{array}{cc}\tan \bar{\theta }x+\frac{\mathrm{qx}\left[2\delta \right(t)-x]}{2\nu \cos \bar{\theta }},\hfill & \mathrm{in}\left[0,\delta \left(t\right)\right],\hfill \\ \mathrm{qt},\hfill & \mathrm{in}\left(\delta \left(t\right),L-\delta \left(t\right)\right),\hfill \\ \tan \bar{\theta }(L-x)\hfill & \hfill \\ +\frac{q(L-x)\left[2\delta \right(t)-L+x]}{2\nu \cos \bar{\theta }},\hfill & \mathrm{in}\left[L-\delta \left(t\right),L\right],\hfill \end{array}\right.$$

with δ(t) given by

$$\delta \left(t\right)=\frac{\sqrt{{\nu }^2{\sin }^2\bar{\theta }+2\nu q^2\cos \bar{\theta }t}-\nu \sin \bar{\theta }}{q},$$

until

$$t\approx \left(\sin \bar{\theta }+\frac{L}{4\nu }\right)\frac{L}{2\cos \bar{\theta }}.$$

Thereafter, the trapezoidal shape takes the almost triangular shape in Equations (15) and (16).

Figure 3.

Evolution of the sand pile profile on a table in the presence of a wall. In (a,b), there is an open end on the left and a wall on the right. In (c,d), there are walls on both sides. Sand is constantly poured with q = 0.05 in the sub-intervals $D_w$ highlighted by the dashed segments in the upper part of the graphs. Specifically: (a) $D_w=$ (0.3, 1); (b) $D_w=$ (0, 0.5); (c) $D_w=$ (0.45,0.55); (d) $D_w=$ (0.25, 0.35) ∪ (0.75, 0.85). In (a, b) the stationary profile is shown in red. In (c, d) there is no stationary profile.

Note that when the input interval is not the entire domain, the peak does not always correspond to the center of the location of the source term as observed in [16]. In fact, focusing for instance on Figure 3(b), the profile starts having a trapezoidal shape with the oblique sides having an angle equal to the angle of repose (fixed here to 45° for easier comparison with the results of the simulations obtained using a two-layer model [26]). One of the sides ends with the end of the open table and it is fixed. This means that after some time the shape becomes an isosceles triangle. The tip gradually moves to the right because part of the sand poured on the top of the triangle moves to the left side and drops down the plate while the part moving on the right side remains on the table. This continues until the tip reaches the end of the interval of the source because at this point all the sand is poured on the left side of the triangle dropping down the plate. Similar dynamics hold for Figure 3(c,d). In all of them, the coefficient ν is normalized to 1 corresponding to a suitable rescaling of time. The simulation set-ups were chosen to compare the results with those obtained by Falcone and Finzi Vita [26] using the two-layer model. In this comparison the reader must be aware that our h should be compared with the sum of the resting and the motile layer. If the former is very similar to our results (especially the minimal solution of the two-layer model), the latter has sometimes unphysical behaviors.

With respect to Figure 3, in Figure 4(a,b) there is a wall at the right extremum of the interval. In Figure 4(a), the profile starts off flat close to the wall having an oblique slope close to the other extremum of the source domain. In Figure 4(b) the evolution is similar to the case in Figure 3(b), but for the fact that when the sand reaches the wall it stops and does not drop down the plate. However, no further sand accumulates in the region [0.5, 1] after the tip of the triangle has reached x = 0.5 because all the sand flows to the left. We note that, as shown in the careful numerical results presented in [25], in two-layer models the height of the sand pile on the right of the input interval continues to grow (see their Figure 1).

Figure 4.

Stationary profile of a sand pile on a square open table due to a constant source q∕ν = 0.5 located in the subdomains: (a) $D_w=$ (0, 1) × (0,1); (b) $D_w=$ (0.45, 0.55) × (0.45, 0.55).

In Figure 4(c,d) there are walls on both sides, which implies that there are no stationary solutions, but the “container” gradually fills. In fact, in the literature, these correspond to silo filling problems in one dimension. We note that for the set-up of Figure 4(c), Finzi Vita [28] noted an unphysical “small depression in the central region corresponding to the source support,” which is not present here. Solutions of the two-layer models, though supported by analytical results, are smoother then expected. Conversely, the model presented here has no problem in generating and handling the well-known formation of corners.

In Figures 5 and 6, two-dimensional simulation results are shown in different domain configurations. In Figure 5(a), the classical pyramidal shape for the open table problem is obtained when sand is poured everywhere, while in Figure 5(b) it is only poured in a central square region, so that only the central portion of the sides of the boundaries is reached by sand (and the slopes flattens up) and not the corners.

Figure 5.

Two-dimensional simulations in the presence of walls and obstacles with constant source q∕ν = 0.5 (a,c), q∕ν = 0.05 (b), and q∕ν = 0.1 (d) in the subdomains: (a, b) $D_w=$ (0, 1) × (0,1); (c) $D_w=$ (0.35, 0.65) × (0.7,0.95) (between the two columns); (d) $D_w=$ (0.4, 0.6) × (0.4,0.6).

Figure 6.

Set up of the experiments (a) initial state; (b) final state.

As in the one-dimensional cases, the results of the simulations can be compared with those obtained using the two-layer models in [23, 25], always recalling that the height of the pile is the sum of the heights of the static and the rolling layers in two-layer models. In particular, the final pile surfaces of Figure 5 can be compared with those in [25], while results in Figure 6(a,b) can be compared with those obtained in [28]. We note that in this case the solution of the two-layer model presents peaks for the rolling layer, and therefore for the total height of the sand pile, at the extrema of the walls where there is a switch between a Dirichlet and a Neumann boundary condition (see their Figures 3–5). Such an undesired effect is not produced by the DPSM proposed here.

Figure 6 (c,d) show the behavior of the DPSM in more complex configurations. We note in particular the formation of discontinuities in the slope on the side of the column opposite to the one where sand is poured. In fact, there the two flows going around the opposite sides of the column(s) meet with a critical slope.

5. Experimental validation: axisymmetric dam-break granular flows

In this section, we validate the mathematical model through a comparison with experimental results from different authors. The experiment performed by Lajeunesse et al. [29], Lube et al. [30], Thompson and Huppert [35], and Roche et al. [36] consists of filling a cylindrical tube of radius $r_0$ with a granular material up to a height $h_0$ to form a cylindrical column (see Figure 7). To trigger the collapse suddenly the tube is raised vertically. The collapse starts from the edges of the sand column moving under free fall developing a front flow at the base of the sand pile that propagates outwards. For sake of completeness, we recall that similar experiments in a 2D-cartesian set-up were performed in [31, 32, 33, 34]. Here sand is accumulated in a channel closed at one hand by a wall and at the other hand by a gate that is impulsively removed. The dynamics is very similar to the axisymmetric case, but for the fact that in this experiments there might be wall effects. For this reason, we selected the axisymmetric set-up. This kind of column/stepcollapse experiments (two-dimensional and axisymmetric case) has been reproduced numerically using Lagrangian approaches, such as Discrete Element Method [38, 39, 40, 41] and Material Point Method [42].

Figure 7.

Comparison of the evolution of dam-break between experiments from Lajeunesse et al. [27] (left column) and simulations obtained on the basis of DPSM (right column): (a) a = 0.56 giving rise to a truncated cone; (b) a = 0.8 giving rise to a cone, (c) a = 5.4 with strong inertial effect not included in the DPSM.

Lajeunesse et al. [29] discussed the results as a function of the aspect ratio of the column, defined as $a=h_0∕r_0$ . In particular, they found that for a < 0.74 the final deposit results in a truncated cone of height close to h₀ and the angle of the sides is close to the angle of repose (that is reported to be close to 21°). For a > 0.74 there is no upper surface and the final configuration is a cone with height smaller than h₀.

The critical value 0.74 is not surprising. In fact, a cylinder with aspect ratio equal to 0.74 has the same volume as a cone with the same height and base angle nearly equal to 23°. The small discrepancy is due to the fact that a free fall of grain and inertial effects are always present, especially at larger aspect ratios and these phenomena are not included in the model. In addition, while in the cylinder the sand packing is near the close-packing value, in the cone the sand might not be so compact.

It is also observed that when the aspect ratio grows (say α > 3) the dynamics related to grain falling become prevalent with respect to grain sliding. Actually, for aspect ratios much larger than one, the entire column falls down giving rise to a Mexican hat shape with angles that are shallower than the angle of repose due to stronger inertial effects and a dynamics that is different from that studied in this paper. For this reason, we actually compare the outcomes of the degenerate parabolic model only with experimental data in the range 0 ≤ a ≤ 2.

A comparison of the experimental results and the outcome of the model proposed here is plotted in Figure 8, which shows the outcome of the simulation for a = 0.56, 0.8, 5.4, together with a reproduction of the experimental snapshots in [29]. From a qualitative point of view, the final configurations obtained by means of computational simulation are very close reproductions of the experimental ones, as well as intermediate non-equilibrium steps in the first two cases. The third one is not because the entire column collapses down and actually there is no creep layer.

Figure 8.

Normalized final height (a) and radial displacement (b) of the pile for different initial aspect ratios. Comparison between the numerical simulation of the degenerate parabolic model and experiments by Lajeunesse et al. [28], Lube et al. [29], and Roche et al. [35].

Lube et al. [30] performed the same experiment using sand, salt, and other different granular materials, which have angles of repose between 30° and 35°. Consequently, the value that separates the truncated cone from the cone regime is different. In spite of the variability of the angle of repose, they give a unique range for all the cases: a < 1 for truncated cone shapes and a > 1 for cone shapes. Roche et al. [36] performed this experiment using coarse and fine granular material as well as dry and fluidized media. In the light of model-validation purposes, we only consider the dry cases. Thompson et al. [35] carried out experiments with a wide range for the aspect ratio a, giving just a few values in our range of interest.

In Figure 9 the outcome of the degenerate parabolic model is compared with the experiments mentioned above focusing on the characteristic dimensions of final pile configuration for different aspect ratios. In particular, Figure 9(a) shows the normalized final height of the pile $h_f$ . With regards to the degenerate parabolic model, 21° and 30° (namely DPSM-21 and DPSM-30) are used on purpose to provide a comparison with the different experiments. Both of them precisely switch from the truncated cone to cone regime (at a ≈ 0.7 for $\bar{\theta }=21^{°}$ , and at a ≈ 1 for $\bar{\theta }=30^{°}$ ). Considering the cone range, both of them capture the $h_f$ -rate of change decrement due to surface lowering.

Figure 9.

Normalized potential energy dissipation for different initial aspect ratios, comparison between numerical simulation, experimental measurements, and semi-empirical law (17).

Figure 9(b) shows normalized radial displacement. The DPSM-21 curve matches accurately its empirical counterpart (specifically, Lajeunesse et al. [29]), even though for high aspect ratio the effect of free fall dynamics beginning to be visible. In fact, as reported by Lajeunesse et al. [29] and Lube et al. [30] in this regime, the final surface slope is smaller than the angle of repose, leading to higher final radius. Surprisingly, the results of Lube et al. [30], dealing with $\bar{\theta }=30^{°}$ , should lead to radial displacement smaller than $\bar{\theta }=21^{°}$ experiments by Lajeunesse et al. [29]. Conversely, they obtained slightly larger radial displacement, while the DPSM-30 curve is coherently below the DPSM-21 trend. A further experiment with different angles of repose could clarify this issue.

Tables 1 and 2 summarize the results and provide fitting curve formulas of experimental and numerical trends. As far as the degenerate parabolic model is concerned, data have been fitted with linear curves, differentiating according to the aspect ratio regimes defined previously.

Table 1.

Summary of normalized $h_f$ trends from experiments and DPSM.

	$\bar{\theta }$	$\frac{h_f}{r_i}$
Lajeunesse et al. [29]	21°	$\begin{array}{c}a,a\le 0.74\hfill \\ 0.74,a>0.74\hfill \end{array}$
Lube et al. [30]	30°	$\begin{array}{c}a,a\le 1\hfill \\ 0.88a^{1∕6},a>1\hfill \end{array}$
Roche et al. [36] - coarse	25°	$\begin{array}{c}a,a\le 0.7\hfill \\ 0.7,a>0.7\hfill \end{array}$
Roche et al. [35] - fine	28°	$\begin{array}{c}0.84a,a\le 1\hfill \\ \approx 0.95,a>1\hfill \end{array}$
Thompson [35]	30°	0.74a
DPSM-21	21°	$\begin{array}{c}a,a\le 0.7\hfill \\ 0.186a,a>0.7\hfill \end{array}$
DPSM-30	30°	$\begin{array}{c}a,a\le 0.95\hfill \\ 0.244a,a>0.95\hfill \end{array}$

Table 2.

Resume of normalized $r_f$ –r₀ trends from experiments and this work: $f_1=1∕\left(2\tan \bar{\theta }\right)(a+\sqrt{4{\tan }^2\bar{\theta }-a^2∕3})-1$ ; $f_2=\sqrt{3a∕0.74}-1$ .

	$\bar{\theta }$	$\frac{r_f-\mathrm{ri}}{r_i}$
Lajeunesse et al. [29]	21°	$\begin{array}{c}{f}_1,a\le 0.74\hfill \\ f_2,a>0.74\hfill \end{array}$
Lube et al. [30]	30°	$\begin{array}{c}1.24a,a\le 1.7\hfill \\ 1.6a^{1∕2},a>1.7\hfill \end{array}$
Roche et al. [35] - coarse	25°	$\begin{array}{c}1.40a,a\le 1.6\hfill \\ 1.82a^{1∕2},a>1.6\hfill \end{array}$
Roche et al. [36] - fine	28°	$\begin{array}{c}1.61a,a\le 1.6\hfill \\ 1.84a^{2∕3},a>1.6\hfill \end{array}$
This work 21°	21°	$\begin{array}{c}1.068a,a\le 0.7\hfill \\ 0.554a,a>0.7\hfill \end{array}$
This work 30°	30°	$\begin{array}{c}0.733a,a\le 0.95\hfill \\ 0.449a,a>0.95\hfill \end{array}$

Lajeunesse et al. [29] also provided potential energy dissipation during the dam-break for different aspect ratio experiments. We compare their empirical trend with the degenerate parabolic model using $\bar{\theta }=21^{°}$ . In Figure 10, numerical results follow in a satisfactory way Lajeunesse et al.’s semi-empirical fit

$$\frac{\Delta E}{E_i}=\left\{\begin{array}{cc}\frac{a}{6\sqrt{3}{\tan }^2\bar{\theta }}\sqrt{12{\tan }^2\bar{\theta }-a^2}\hfill & \mathrm{for}a<0.74;\hfill \\ 1-\frac{0.74}{2a}\hfill & \mathrm{for}a>0.74.\hfill \end{array}\right.$$

(17)

Figure 10.

Application of the DPSM to dune dynamics.

Again, for high aspect ratios, free-fall dynamics lead to a slight difference in the final configuration and in the dissipated energy accordingly.

6. Discussion

In this work, we have deduced a degenerate parabolic model for sand slides based on continuum mechanics methods. It is assumed that the thickness of the creep layer is small, that the sand density in it is close to close packing, and that grains move in the direction of the steepest descent with the velocity limited by friction and a drag force that mimics the interaction among grains. Owing to the above assumptions, the model cannot describe the massive motion of granular materials in their generality, such as the one recently studied by Edwards and Gray [43].

However, despite its simplicity, the model is able to catch all the classical behaviors of sand piles in different conditions, such as starting from supercritical conditions, growing because of external material pouring, in the presence of walls and obstacles or with open boundaries. Each of these particular cases does not require any ad-hoc modifications of the model. From a numerical point of view, different examples have been reported, in one and two dimensions, with different boundary conditions and internal obstacles, in order to show the high flexibility of the model, which actually does not present any problems in addressing these kinds of issues. Moreover, numerical validation by means of comparison with different experimental works on column collapse has been provided. It shows the capability of the model to reproduce experimental results in a satisfactory way, from both the qualitative and the quantitative point of view. if the aspect ratio of the column is not too large because in that case sand grains fall rather than slide and intertial effects not included in the model become important.

It is worth pointing out that the DPSM can be used in dune dynamics applications and sand-bed evolution, where avalanching locally acts to modify slightly supercritical configurations formed by the continuous transport of sand from the upwind to the downwind side of the dune. Referring to Preziosi et al. [1] for more details on how the fluid dynamic coupling could be performed, in Figure 11, we show an example of how a transverse dune with a slightly supercritical slip face would evolve to reach a stable configuration.

All the previous characteristics make the degenerate parabolic model very interesting to describe sand avalanches, in particular when coupling between multiphase turbulent fluid dynamics and the dynamics occurring at the sand surface is required. In fact, the model meets the requirements of describing the evolution of the sand surface in an accurate way, but at the same time it is computationally fast, which is a relevant plus considering the highly demanding multiphase simulations.