English (pdf)
Article in xml format
Article references
Automatic translation
Send this article by e-mail
Cited by SciELO 
Cited by Google 
Similars in
SciELO 
Similars in Google 
Permalink
Introduction
The study of granular avalanches provides a basis for understanding a variety of mass movements, such as rock avalanches, snow avalanches, debris flows and pyroclastic flows Denlinger and Iverson (2001). Some authors (Denlinger and Iverson, 2001, Savage and Hutter, 1989, Gray et al., 1999, Wang et al., 2004, Pastor et al., 2015, Koch et al., 1994) have made efforts to perform laboratory experiments and numerical simulations of granular avalanches. The continuum mechanics framework supplies a helpful theoretical standpoint to represent multiple types of flows by using the laws of conservation of mass and momentum. A reduced form to use the conservation laws applied to granular flows is by means of a shallow water approach. The bidimensional equations that describe the movement of a granular mass taking into account the stresses, caused by either the grains or the fluid.
The granular flow equations are prone to develop discontinuities due to its nonlinear nature. Often, the solution to these equations can be obtained by mean of numerical methods such as finite differences, finite volumes, and finite elements, known as low-order techniques. Results using methods of low order are reported by Gray et al., (1999)), Hutter and Koch, (1991, Hutter and Greve, (1993), Hutter et al., (2005), Ouyang et al., (2013), Savage and Hutter, (1989), Tai et al., (2002), Wieland et al., (1999), Denlinger and Iverson, (2001), Denlinger and Iverson, (2004), Koschdon and Schafer, (2003), Vollmoller, (2004), Wang et al., (2004)), where the numerical solution shows a similar structure when compared with experimental data. However, the techniques mentioned above present significant numerical error, such as artificial diffusion and spurious oscillations that affect the stability, accuracy, and conservative properties of the solution as is shown by Fletcher, (1991). In fact, Diamessis and Redekopp, (2006) mentioned that numerical stability resulting from low-order finite-difference schemes is caused by the well-known inherent artificial dissipation, which might generate a significant error in the longtime integration. This issue has been handled with a series of robust methods known as shock-capturing schemes such as ENO, WENO, NOC, TVD, HLLC, among others (Denlinger and Iverson, 2001, Tai et al., 2002, Wang et al., 2004, Toro, 2013). However, it is still necessary to perform significant efforts for developing numerical techniques that allow solving the systems of equations as well as possible to describe these natural phenomena (Castro-Orgaz et al., 2015, Iverson, 2014).
In order to minimize the error associated with numerical diffusion generated by using low-order techniques, researchers as Hesthaven and Gottlieb, (1996), Reed and Hill, (1973), Diamessis and Redekopp, (2006), Giraldo and Warburton, (2008) developed element-based high-order discontinuous techniques. Benchmark cases with an exact solution were used to evaluate the accuracy and convergence of spectral methods by solving the inviscid shallow water equations (Escobar-Vargas et al., 2012). Based on the evidence from the high-order precision of these methods 0(10-11), the potential of an element-based spectral method for solving a set of equations that describe the movement and deposition of granular avalanches is evaluated, taking into account its lack of numerical diffusion. This work was performed with the aim of test the potential of high-order methods for more applied cases since these techniques have been implemented in idealized cases.
The numerical method implemented is based on a collocation approach of multiple non-overlapping quadrilateral subdomains, which are connected by a penalty term that ensures the stability of the solution by imposing a weak continuity at the subdomain interfaces (see Appendix A). This technique is called spectral multidomain penalty method, henceforth SMPM (Hesthaven and Gottlieb, 1996, Escobar-Vargas et al., 2012). In order to maintain the stability of the spatial discretization technique, a strong-stability preserving Runge-Kutta was implemented for the temporal derivative (Appendix B), simultaneously, a spectral filter was imposed to control the numerical oscillations (Appendix C). The potentiality of the SMPM is evaluated in two ways. The first one is by means of a grid independence test that shows the error in mass conservation introduced by the numerical scheme with the spatial refinement. The grid independence test shows that the error in mass conservation does not improve as much as is reported in the literature for the case of a granular avalanche. The second means of evaluating the SMPM involves comparing the numerical simulations and the experimental measurements performed by Denlinger and Iverson, (2001). The comparison shows that the SMPM provides comparable results for depths of the flow concerning the experimental measures and other numerical techniques. Finally, a discussion of boundary conditions, lateral earth pressure coefficients and the momentum correction factor were raised due to its incidence when solving the granular flow equations.
The paper is organized as follows. In section 2, the granular flow equations are presented. In section 3, the numerical setup for the case of study, such as initial and boundary conditions, and treatments are described for zones with or without material. The analysis of a grid independence test is performed in section 4; and in section 5, the comparison between experimental and numerical simulation is performed where boundary conditions, lateral earth pressure coefficients and momentum correction factor are analyzed. A discussion and conclusions are presented in section 6 and 7, respectively. In Appendices A, B, and C we present the numerical techniques with which the equations are discretized and stabilized.
Governing equations of the granular flow
Many granular avalanche models have been proposed based on the 2D shallow water approach as shown by Savage and Hutter, (1989), Iverson and Denlinger, (2001), Denlinger and Iverson, (2004), Hutter et al., (2005), Wang et al., (2004), where the main differences are focused on evaluating some features of the stresses. In this sense, the model proposed by Iverson and Denlinger, (2001) was selected as a sample in order to evaluate the potential of the SMPM for solving granular flow equations (see Appendix A). The temporal derivative is handled by a strong-stability preserving Runge-Kutta (See Appendix B). The model can be presented in conservative form as
where q denotes the conservative variables vector. F(q) and G(q) represents the horizontal fluxes in the x and y direction, respectively, and S(q) is the vector that contains the source terms, which dissipate the energy (Escobar-Vargas et al., 2012). By decomposing this system of equations we have
In this set of equations, h is the flow depth, m
x
= hū and my =
, where ū and
are the averaged velocities in x and y directions. g
x
= g sin θ, g y = 0 and
g
z
= g cos θ are the gravitational acceleration components in x, y and z directions, respectively, g = 9.81m / s is the gravitational constant and θ is the inclination angle; λ = P
bet
/ ρg
z
h is the ratio of the basal pore fluid pressure to the total normal basal stress. λ is evaluated from an analytical summation expression that is a function of hydraulic diffusivity D and the initial pore pressure ratio k such as is shown in Iverson and Denlinger, (2001). k
a
/
p
is the lateral earth pressure coefficient; and β is the momentum correction factor that comes from the integration of the advective terms of the Navier-Stokes equations (Zhou, 2004). The structure of equations (2-5) is similar to the one of inviscid shallow water equations, where the source terms S
x
and S
y
take into account the stresses generated by fluid and solids. Details of these features are shown in equations 6 and 7.
where Ф int is the internal friction angle of displaced materials, Ф bed is the friction angle of bed material, n f is the fluid volume fraction (i.e., the porosity), µ is dynamic viscosity of the fluid, and θ x , θ y are the components of local angle on the bed in the x and y directions, which are described in Section 3.1. The lateral earth pressure coefficient (k a / p ) that appears on Equations (3-7) can be represented by
This coefficient associates the horizontal and vertical stresses (Savage and Hutter, 1989; Iverson and Denlinger, 2001). Equation 8 shows that k a / p has a discontinuous nature, indicated by the plus/minus sign, which jumps from active to passive states (i.e., ∂u / ∂x > 0 or ∂u / ∂x < 0 for the x direction and ∂v / ∂y > 0 or ∂v / ∂y < 0 for the y direction). This feature makes the solution to be prone to become unstable. A better formulation of this coefficient is given by a regularized version of equation 8 [Tai and Gray, 1998]. The discontinuity is avoided by using the following expression:
where A is assumed to be a monotonically decreasing functions (Equation 10), which depends of the gradient of the velocity.
In Equation 10, ψ determines the steepness of the transition between active and passive coefficient, and ∂u / ∂x0 is chosen as the center of the transition (i.e., when ∂u / ∂x = 0 ) (Tai and Gray, 1998). This regularized lateral earth pressure coefficient is applied in both the x and y directions and yields a more stable solution (see section 5.2). Another version of the lateral earth pressure coefficient was implemented in Tai et al., (1999), Gray et al., (1999),Tai et al., (2001), Wieland et al., (1999), Gray, (2003), where the relationship between vertical and horizontal stresses in the y direction is a function of the lateral earth pressure coefficient in the x direction.
In addition, a constant lateral earth pressure coefficient (k nat ) was employed assuming that internal and bed friction angles are equal (φ int = φ bed ). The main importance of implementing this constant coefficient is to perform the evaluation of other parameters in a more controlled manner by avoiding the discontinuous nature of equation 8 or even the softened function 10, by using equation 11 (Iverson and Denlinger, 2001).
The set of equations 1 - 5 belongs to a family of models in which one set of equations is employed, and the stresses are split into two components, solid and fluid. This characteristic makes a significant difference since it avoids the difficulties from a real two-phase model, where a set of equations is needed for each phase and its coupling (Pudasaini, 2012). The fraction of the fluid is a constant parameter, meaning that the proportion of fluid is homogeneous, and it does not vary in time. The role of the fluid fraction in the configuration used here (i.e., granular ow) seems to be irrelevant since the dynamic viscosity of the air is a small value (i.e., 2x10-5). Hence, stresses produced by the fluid are weak regarding the stresses caused by the solids.
Geometrical basis and initial conditions of the granular flow
Denlinger and Iverson (2001) presented a small-scale experiment of a dry granular avalanche. They built a 1.2 m long rectangular flume with an inclined bed at 31.4° with respect to the horizontal surface. There is a curved section with a 10 cm radius of curvature between the inclined and the horizontal surface (Figure 1). Dry quartz sand with a grain diameter of 0.5 mm was used. The initial setup of the sand is wedge-shaped and is located behind a vertical gate positioned at a 62.5 cm slope down from the top of the flume, which is discharged suddenly when the gate opens (Figure 1a) (Denlinger and Iverson, 2001). The flume and gate width is 20 cm. The geometry of the bed varies only in the x direction. This variation can be written as
Figure 1 Initial conditions employed in this work. a) 3D view of the topographic basis of an inclined and curvilinear plane, with a wedge-shaped initial condition. Additional initial conditions, b) cylinder and c) Gaussian bump of sand.
Based on Equation 12, the geometrical configuration of the flume is split into three parts: the slope (i.e., 0 ≤ x ≤ x l ), the curved zone (i.e.,), and the horizontal zone (i.e, x¡ ≤ x ≤ x r ) (Figure 1a). In Equation 12, b(x) describes the elevation of the flume, b , is the maximum height of the curved zone, at left side of the curved part, x t is any point of x, on the curved zone (rx), x r is any point of x on rx, and m 0 is the slope of the inclined surface.
It is possible to obtain the slope angle from the geometrical configuration of the basis (i.e., θ x = tan-1 (∂b/ ∂x))and the curvature (i.e., ∂ x / ∂x) that appears in equations 6 and 7. The initial condition was built as a wedge (Figure 1), which is described by the following equation:
The initial condition (h i ) is split into three points, left (x l ), central (x c ) and right point (x r ). The wedge that represents the initial condition is constructed between the left and right points. m 1 and m 2 are the slope of each side of the triangle.
Table 1 shows the properties of the materials used in the numerical simulation presented in this work (Denlinger and Iverson, 2001).
Table 1 Value of the properties of the materials employed in the numerical modeling. Taken from (Denlinger and Iverson, 2001).
| Parameter | Unit | Value |
|---|---|---|
| Basal friction angle ( Ф bed ) | degrees | 29 ± 1.4 |
| Internal friction angle (Ф int ) | degrees | 40 ± 1.0 |
| Solid volume fraction (n s ) | 0.6 | |
| Fluid volume fraction (n f ) | 0.4 | |
| Fluid viscosity (µ) | Pa-s | 2x10-5 |
| Solid density (ρ s ) | kg/m3 | 2650 |
| Fluid density (ρ f ) | kg/m3 | 1 (air) |
| Mixture bulk density (ρ) | kg/m3 | 1600 |
| Hydraulic permeability (D) | m2/s | 0.05 |
| Initial pore pressure ratio (k) | 0.0 |
Boundary conditions
An appropriate implementation of boundary conditions is one of the most important aspects to obtain a robust numerical approximation of any partial differential equation. The relevance of boundary conditions is presented in Section 5.1, where it is shown that different solutions are obtained by using different treatments on boundaries. In order to analyze the influence of the boundary conditions in the solution, two types of boundary conditions were implemented and are presented below.
Figure 2 shows the representation of two boundary conditions. Reflective boundary conditions are imposed in both cases, on boundaries perpendicular to the flow direction (i.e., the velocity vector is
ūn =
0, where
ū = (ū
,
)). In this case, the right boundary must be far enough from the mass movement to avoid any influence in the dynamic of the numerical experiment. Thus the mass can stop by the frictional force solely and not by the effect of an obstacle. Also, slip boundary conditions were imposed on boundaries parallel to the flow direction,
ūn =
0 (Figure 2a). In the second case, "non-slip" boundary conditions were implemented on boundaries parallel to the flow direction, which consists of imposing a velocity value equal to zero,
ū =
0 (Figure2b).
The tracking of the free surface
To solve the granular avalanche equations (Equation 1), it is necessary to define a method to treat regions without material (i.e., h = b). Shallow water equations are defined by a height such that h = H
o
+ ξ, where H
0
is an average depth and ξ is the displacement of the free surface (Kundu and Cohen, 2008). From solving the system of equations, the velocities ū and are computed by means of conservative variables with Equation 14.
The use of average depth H
0
allows the appearance of undesired oscillations below the average height of the flow that cannot appear in dry zones (h = b), when a rigid bed exists. In this work, the average depth is H
0
= 0 for the granular flow equations, then h = ξ. In areas without material (i.e., h = 0) the solution of the set of equations in the conservative form generates a mathematical indetermination (i.e., division by zero to compute ū and
, see Equation 14). In order to avoid this problem, a tolerance value was defined (e.g., Г = 10-4). If h < Г, a trivial solution is imposed (i.e., ū =
= 0) (Bradford and Sanders, 2002; Pudasaini and Hutter, 2007). Otherwise, if h >Г, the velocities are computed without imposing any restriction (Pudasaini and Hutter, 2007).
Stability condition and mass conservation
The stability condition for the SMPM is given by the Courant number, which is defined by Giraldo and Warburton (2008) as follow:
where At is time step size,
is the magnitude of average velocities into the cell, gh is average geopotential height into each cell and
is the grid spacing. The sum of
is known as the characteristic velocity at which the wave travels (LeVeque, 2002). For SMPM, the maximum number of Courant that ensures the stability is CFL ≤1 (Escobar-Vargas et al., 2012).
Conservation of mass was used as the main criterion by which to evaluate the conservative properties of the method. The error introduced by the numerical method was calculated with respect to the initial mass. The mass of the system can be calculated as
where Q represents the domain. The metric to evaluate mass conservation is based on the relative error, defined with respect to the initial values of mass as
In equation 17, RM is the relative error of the mass, and M0, Mt are the corresponding values for the mass at the initial time (t = 0) and final time (t) of the simulation. Thus, the loss or gain of mass is specified at each time step. The stopping criterion was focused on the experimental data since the last measurement available to compare is at 1.5 seconds. In addition, the variation of the position of the centroid of the mass and its velocity were kept below 0(10-2) as a secondary criterion. By another hand, the error in energy conservation was not verified in this work because the source terms in the momentum equations are responsible for the energy dissipation.
Grid independence test
The grid independence test was performed to determine the optimal degree of grid refinement at which the solution does not change when the mesh is reduced. The error in mass conservation (equation 17) was chosen as a reference quantity to assess the differences of the solutions among multiple grids. The mesh refinement was handled with two different strategies. The first option is developed by dividing the domain into multiple subdomains; this is called h-refinement. The second option keeps the numbers of subdomains fixed while increasing the polynomial degree; this is called p-refinement (Boyd, 2001).
Figure 3 shows the results of error in mass conservation for the test case with the two types of boundary conditions (slip and non-slip) described in Section 3.2. This was performed in order to observe whether the error is increased with respect to the slip boundary conditions when the non-slip boundary conditions are imposed. The simulation was performed for 9 cases of h-refinement, where the subdomains were kept square (Figure 3), and 19 cases of p-refinement (i.e., 2 ≤ N ≤ 20). It is possible to observe that the error in mass conservation does not improve when the refinement of the mesh is increased. It is noticeable that the non-monotonicity of the grid independence test disappears when p-refinement is greater than 16th degree, independently of the h-refinement for both boundary conditions. The non-monotonic behavior presented in the grid independence test makes more difficult the choice of an optimal spatial discretization. In addition, the error does not significantly decrease, as is reported with benchmark cases (Escobar-Vargas et al., 2012).
Figure 3 Mass conservation error with a) slip and b) non-slip boundary conditions. (Sub x and Sub y represents the subdomain numbers in x and y directions, respectively).
In addition, the implementation of the collapse of a cylinder and a Gaussian bump of sand for two different average depths (H = 1 and H 0 = 0) (Figure 1b and 1c, respectively) were performed. Considering the Gaussian bump (Figure 4a), it was obtained a very low numerical error 0(10-9) both with and without average depth, whereas the cylinder generates a greater numerical error (Figure 4b). The tracking of the free surface provides an additional error caused by the zero-average depth (H 0 = 0). Otherwise, an average depth greater than zero (H = 1) helps to give a lower error. If the value of the tolerance for the "wet/dry" condition is 10-4, the numbers after the fifth decimal place are neglected when the velocities are computed. These velocities are inputs for the next time-step, which affects the accuracy. Although it is natural that hyperbolic equations develop sharp fronts, with these set up cases it was possible to observe that discontinuities might not mature. This is because of the physics of the modeled phenomena included by means of the source terms since these can be stopped soon.
Boundary conditions analysis
Figure 5 presents the comparison between the experimental data (D-I2001) (Denlinger and Iverson, 2001) and the numerical prediction given by the SMPM with a grid that showed the lower error in mass conservation obtained from the grid independence test (Figure 3). Also, slip boundary conditions and the constant lateral earth pressure coefficient were imposed.
Figure 5 a) Granular flow experiment. b) Numerical prediction given by SMPM for spatial refinement with low numerical error (N = 8, Sub x = 108 and Sub y = 18), slip boundary conditions, and constant lateral earth pressure coefficient.
The numerical prediction is akin to the experiment, particularly in the deposit of the granular flow. The maximum heights of the deposit given by the numerical prediction are approximated to the maximum heights of the experiment with differences of 2 mm, approximately. However, the movement of the material is not well captured, because the "adhered" material on the lateral walls does not appear in the numerical prediction as was observed in the experiment.
Figure 6 presents the comparison between the experimental data (D-I2001) and the numerical prediction given by the SMPM with non-slip boundary conditions and the constant lateral earth pressure coefficient. Non-slip boundary conditions allow us to represent the tail of the granular flow while the granular mass is going downhill. The "adhesion" of the material on the lateral walls given by the numerical prediction is similar to the experimental data (the tail of the flow). However, a deposit zone is generated caused by these boundary conditions (the tail of the flow is still present), which is not the case in the experiment.
Figure 6 a) Granular flow experiment. b) Numerical prediction given by SMPM for spatial refinement with low numerical error (N = 8, Sub x = 108 and Sub y = 18), non-slip boundary conditions, and constant lateral earth pressure coefficient.
Denlinger and Iverson, (2001) solved this set of equations by means of a method called Harten-Lax-van Leer-Contact (HLLC) (Toro, 2013). A quantitative comparison among the numerical methods with respect to the experimental data is presented for five different stages (Table 2). The SMPM allows us to get a comparable solution with the HLLC method, since the error in mass conservation has the same order of magnitude for the case of study, with respect to the experimental data.
Table 2 Euclidean norm of the error for slip and nonslip SMPM, and for HLLC, regarding the experimental data.
| Time (s) | SMPMslip | SMPMnon-slip | HLLC |
|---|---|---|---|
| 0.01 | 0.2227 | 0.2197 | 0.4428 |
| 0.32 | 0.1634 | 0.1583 | 0.1796 |
| 0.53 | 0.1694 | 0.1733 | 0.2108 |
| 0.93 | 0.5699 | 0.6237 | 0.4904 |
| 1.50 | 0.4761 | 0.5739 | 0.4157 |
It is possible to see from the Table 2 that the SMPM has generated less differences with respect to the experiment in the first three stages (i.e., while the material is descending the inclined plane), whereas the HLLC method produced better results in the last two stages (i.e., when the material is forming the deposit, and finally is stopped).
Another comparison between numerical results and experimental data is presented in Figure 7. The continues and discontinuous lines represent two ways to compute the centroid of the mass and the centroid of the velocity field, such as is shown in Pudasaini and Hutter, (2007, p. 149). The asterisk represents the centroid of the experiment computed from the cloud of points as
Figure 7 Comparison between numerical results and experimental data of the centroid of the moving mass and the velocity field. a) The position of the centroid of the descending mass at each time step. b) Centroid of the velocity field. V is referring to the total volume.
It is possible to observe that the numerical prediction matched in a good way the position of the center of the moving mass computed from the experimental data, however, the centroid of the velocity field shows difficulties that the prediction fits the experiment.
Analysis of lateral earth pressure coefficient (k a / p )
The lateral earth pressure coefficient ( k (x y)a/p ) relates the vertical stresses with the horizontal ones of the granular flow. Figures 5bb and 6b show the solution generated with a constant k nat (equation 11), whereas Figure 8b presents the solution given by the slip boundary conditions and a regularized k (x y)a/p (equations 9-10). This specific coefficient tends to produce asymmetric solutions caused by nonphysical oscillations, and hence, steep localized gradients from the velocity can be produced. Since the derivative of the velocity is an input of the lateral earth pressure coefficient, this has repercussions in the behavior of this coefficient and to the next time step. These oscillations can be mitigated either by decreasing the order of the polynomial or by handling the spectral filter (Appendix C).
Fig. 8 a. Granular flow experiment. b. Numerical prediction given by SMPM when the spatial refinement presents low numerical error (N = 8, Sub x = 108 y Sub y = 18), slip boundary conditions and regularized lateral earth pressure coefficient k (x y)a/p .
The use of a regularized k (x y)a/p allows us to represent in a better approach the physics of the processes that occur in these types of flows, by means of the smooth connection between the two states of stress (dilation or compression of the material). By another hand, the use of a constant k nat indicates that neither dilation nor compression zones are presented, but this behavior does not occur in nature. For example, experimental data presented by Pudasaini and Hutter, (2007) showed that velocities can change from one point to another. However, a constant k na was used because it makes it possible to obtain a more stable solution, whereas with a regularized k (x y)a/p the solution is more likely to be unstable, even more with the discontinuous coefficient k (x y)a/p
Analysis of momentum correction factor (β)
The momentum correction factor β (see equations 3 and 4) is a parameter derived from depth-integrating the Navier-Stokes equations, undertaken to obtain the shallow water equations (Zhou, 2004, Trujillo-Vela, 2015). The physical meaning of β is that it modifies the profile of the velocity of the flow. For most of the cases this factor is taken as the unity by assuming that the velocity profile is almost uniform through the depth of the avalanche (Denlinger and Iverson, 2001, Savage and Hutter, 1989). However, this assumption may not be the most appropriate after of the release, in the proximity of obstacles and near to the deposit where both height and vertical velocity may change considerably and, in some cases, abruptly (Pudasaini and Hutter, 2007). For this reason, the effect of the momentum correction factor is analyzed in this section.
Figures 9 to 11 show a longitudinal section at the center of the flume for the experiment (in each Figures a and c) and the numerical results (in each Figures b and d) at five different times (0.10, 0.32, 0.53, 0.93, and 1.50 seconds). Heights are shown on the ordinate and the longitudinal length is shown on the abscissa. The numerical results were computed for five values of β (1.0, 1.05, 1.10, 1.15, 1.20).
Figure 9 x versus h plots for different β at times t =0.10 s and t =0.32 s, and the comparison with the experimental measurements for the same time.
Figure 10 x versus h plots for different β at times t = 0.53 s y t = 0.93 s, and the comparison with the experimental measurements for the same time.
Figure 11 x versus h plots for different β at time t = 1.50 s, and the comparison with the experimental measurements for the same time.
Figures 9a and 9c show the profile of the free surface at the center of the flume of the experimental results at t = 0.10 s and t = 0.32 s, respectively. Figures 9b and 9d show numerical predictions at the same times. The best prediction of the front at t = 0.10 s is given by 3 = 1.0, since the numerical results represent the front as steep as in the experiment. Results with 3 > 1.0 show very blunt fronts (Figure 9b). If 3 = 10.5 at t ≥ 0.32 s, the flow presents more blunt fronts, and the arrival of material at the zone of deposit is sooner, as is shown in Figure 9d.
Figures 10a and 10c show experimental results of the profile of the free surface at the center of the flume at t = 0.53 s and t = 0.93 s, respectively. The best predictions at t = 0.53 s is given by β= 1.15, because this value of β presents a distribution of material similar to the distribution in the experiment. At the same time, if β = 1.0, a significative amount of material has descended the slope and the generation of the deposit begins prematurely, with respect to the experiment. At t = 0.93 s as shown in Figure 10c. The experimental measurement shows that the deposit is almost fully formed (Figure 10c). A good approximation for the tail of the flow and the deposit is reached when β = 1.10 -1.15 (Figure 10b and d), although the tail is longer than in the experiment. Also, Hutter et al., (2005) show this behavior. On the other hand, there is a numerical oscillation on the left side of the deposit when β = 10 _ 1.10caused by the slope change in this zone and the selected β.
Figure 11b shows that the best prediction of the deposit at t = 1.50, which is obtained by β = 1.0 because the tail is less noticeable (i.e., thinner) unlike if other values for the momentum correction factor are used, and the experiment does not show a tail. In addition, the maximum heights of the deposit were predicted better when β = 1.0.
The analysis of the momentum correction factor offers information about changes in the shape of the granular flow prediction by imposing different values of this factor. Therefore, it is possible to obtain a better prediction in the process where the modeler is most interested in, either during the movement or the deposit, by using β > 1.0 or β = 1.0, respectively. Hence, we propose that future research examine this parameter in order to generate a momentum correction factor that can change depending on what period of movement is being examined when the shallow water approach is used.
Discussion
The results obtained in this work allow us to emphasize some important aspects regarding the solution of granular flow equations. The first aspect is the importance of performing a grid independence test for any numerical method employed in order to obtain an optimal grid. The grid independence test is important because the order of the error can vary significantly for each numerical method and test case. Therefore, this test helps us to find refinement of the grid that presents low error in the mass conservation. Another aspect is the importance of specifying details of the numerical setup, such as initial and boundary conditions, stabilization techniques (Appendix C), and "wet-dry" treatments, since most of the times a part of this information is omitted.
Obviating the error caused by the difficulty of solving hyperbolic systems, the following most important issue that causes the error is the tracking of the free surface for zones without material (i.e., wet-dry cells for shallow water equations), as was shown in Section 3.3. None of the boundary conditions used in this work present a completely satisfactory prediction about the distribution of material near lateral boundaries, since the slip boundary condition solely allow us to represent adequately the deposit of material and the non-slip conditions allow us to predict correctly the movement of the mass. This highlights the importance of exploring different types of boundary conditions, with which a great representation of the movement and deposit of the mass can be obtained.
As shown in Section 2, the use of a discontinuous lateral earth pressure coefficient ( k (x,y)a/p ), such as was proposed by Savage and Hutter, (1989), makes prone to get unstable solutions. Furthermore, as shown in Section 5.2, the use of a constant lateral earth pressure coefficient promotes the getting of more stable solutions but is not the best choice from the physical standpoint, because the material suffers thinning and thickening during the movement. Therefore, the regularization of k (x,y)a/p as proposed by Tai and Gray, (1998), is required.
Finally, the tests performed using the momentum correction factor (β) allow us to conclude that the most suitable value of the momentum correction factor is β = 1.0 after release and close to the deposit. On the other hand, to obtain a better representation of the distribution of material during movement, to mitigate the steep slope of the propagation front and to avoid the premature generation of the deposit it is necessary to use β > 1.0. This parameter shows the limitation of the depth-integrated models to represent a wave propagation front when the unity of the 3 is selected, as is commonly performed.
Conclusions
The solution of granular flow equations with a spectral method such as the SMPM, have not been reported yet. Based on the results reported for the first time in this work, the SMPM solves appropriately the granular flow equations proposed by Iverson and Denlinger (2001) and it provides comparable predictions with the experimental data and others' numerical.
We observed from the results that the error in mass conservation increases significantly when a discontinuity is imposed as an initial condition or the sharp fronts are developed during the movement of the mass and, also when the treatment for the tracking of the free surface is used. By another hand, it is possible to use the advantageous property of hp-refinement of the SMPM, instead of h-refinement, as occurs with low-order techniques.
Slip and non-slip boundary conditions modify the mass distribution considerably. Numerical solutions with slip boundary conditions provide better predictions about the deposit, whereas non-slip boundary conditions provide better predictions about the movement of the granular flow.
Instabilities in the solution of the granular avalanche equations introduced by the discontinuous lateral earth pressure coefficient proposed by Savage and Hutter (1989) can be mitigated by the regularization of this coefficient, such as shown in Tai and Gray (1998). Also, it is possible to obtain a stable solution with a constant lateral earth pressure coefficient.
The best prediction of the beginning of the flow and the deposit are obtained when the momentum correction factor is imposed equal to the unity (β = 1.0), whereas, the best prediction of the movement of the flow is obtained if the momentum correction factor is greater than the unity (β > 1.0).
