SciELO - Scientific Electronic Library Online

 
vol.23 issue4Earthquake predictions and scientific forecast: dangers and opportunities for a technical and anthropological perspectiveInversion Technique of Physical Parameters Based on Regularization Extension Depth Constraint author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand

Journal

Article

Indicators

Related links

  • On index processCited by Google
  • Have no similar articlesSimilars in SciELO
  • On index processSimilars in Google

Share


Earth Sciences Research Journal

Print version ISSN 1794-6190

Earth Sci. Res. J. vol.23 no.4 Bogotá Oct./Dec. 2019  Epub Apr 20, 2020

https://doi.org/10.15446/esrj.v23n4.77683 

Originals Articles

A spectral multidomain penalty method solver for the numerical simulation of granular avalanches

Método penalizado de multidominios espectrales para la simulación numérica de avalanchas granulares

Mario Germán Trujillo-Vela1  * 

Jorge Alberto Escobar-Vargas1 

Alfonso Mariano Ramos-Cañón2 

1Departamento de Ingeniería Civil, Pontificia Universidad Javeriana, Bogotá, Colombia.

2Departamento de Ingeniería Civil y Agrícola, Universidad Nacional de Colombia, Bogotá, Colombia.


ABSTRACT

This work presents a high-order element-based numerical simulation of an experimental granular avalanche, in order to assess the potential of these spectral techniques to handle geophysical conservation laws. The spatial discretization of these equations was developed via the spectral multidomain penalty method (SMPM). The temporal terms were discretized using a strong-stability preserving Runge-Kutta method. Stability of the numerical scheme is ensured with the use of a spectral filter and a constant or regularized lateral earth pressure coefficient. The test case is a granular avalanche that is generated in a small-scale rectangular flume with a topographical gradient. A grid independence test was performed to clarify the order of the error in the mass conservation produced by the treatments here implemented. The numerical predictions of the granular avalanches are compared with experimental measurements performed by Denlinger and Iverson (2001). Furthermore, the boundary conditions and parameters such as lateral earth pressure coefficients and the momentum correction factor were analyzed in order to observe the incidence of these features when solving the granular flow equations. This work identifies the benefits and weaknesses of the SMPM to solve this set of equations, and thus, it is possible to conclude that the SMPM provides an appropriate solution to the granular flow equations proposed by Iverson and Denlinger (2001) and comparable predictions for the experimental data.

Keywords: Granular avalanches; Numerical simulation; Spectral multidomain penalty method; Parameter analysis

RESUMEN

Este trabajo presenta la simulación numérica de una avalancha granular a escala experimental basada en elementos de alto orden, con el fin de evaluar el potencial de estas técnicas espectrales para resolver las leyes de conservación en geofísica. La discretización espacial de estas ecuaciones se desarrolló a través del métodos penalizado de multidominios espectrales (SMPM). Los términos temporales se discretizaron utilizando un método Runge-Kutta de preservación de estabilidad fuerte. La estabilidad del esquema numérico se garantiza con el uso de un filtro espectral y un coeficiente de presión lateral de tierras constante o regularizado. El caso de prueba es una avalancha granular que se genera en un canal rectangular de pequeña escala con gradiente topográfico. Se realizó una prueba de independencia de malla para aclarar el orden del error en la conservación de la masa producida por los tratamientos previamente mencionados. Las predicciones numéricas se comparan con las mediciones experimentales realizadas por Denlinger & Iverson (2001). Además, se analizaron las condiciones de contorno, y los parámetros como el coeficiente de presión lateral de tierras y el factor de corrección de momentum para observar su incidencia al resolver las ecuaciones del flujo granular. Este trabajo identifica las bondades y debilidades del SMPM para resolver este conjunto de ecuaciones y es posible concluir que el SMPM proporciona una solución adecuada de las ecuaciones propuestas por Iverson & Denlinger (2001). Además, las predicciones son comparables a los datos experimentales y resultados numéricos dados por otros esquemas.

Palabras clave: Avalancha granular; simulación numérica; método penalizado de multidominios espectrales; análisis de parámetros

Nomenclature 

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.

Granular flow case

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, xx 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).

Figure 2 Reflective boundary conditions, (a) slip and (b) non-slip in the computational domain. 

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.

Figure 4 Mass conservation error in function of time with and without average depth and without a topographic gradient, a) Gaussian bump and b) cylinder. 

Comparison with experimental data

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).

Acknowledgements

The authors acknowledge Richard Iverson and Roger Denlinger for providing valuable experimental and numerical data for comparison with our numerical results. Fundación CeiBA supported this work within the scholarship of doctoral studies of the first author.

References

Blackburn, H. M., & Schmidt, S. (2003). Spectral element filtering techniques for large eddy simulation with dynamic estimation. Journal of Computational Physics, 186(2), 610-629. [ Links ]

Boyd, J. P. (2001). Chebyshev and Fourier spectral methods. Courier Dover Publications, Mineola, New York, 2nd edition, 670 pp. [ Links ]

Bradford, S. F. & Sanders, B. F. (2002). Finite-volume model for shallow-water flooding of arbitrary topography. Journal of Hydraulic Engineering, 128(3), 289-298. [ Links ]

Canuto, C., Hussaini, M. Y., Quarteroni, A., & Thomas Jr, A. (2012). Spectral methods in fluid dynamics. Springer Science & Business Media. Berlin, Germany, 566 pp. [ Links ]

Castro-Orgaz, O., Hutter, K., Giraldez, J. V., & Hager, W. H. (2015). Nonhydrostatic granular flow over 3-D terrain: New boussinesq-type gravity waves? Journal of Geophysical Research: Earth Surface, 120(1), 1-28. [ Links ]

Costa, B. & Don, W. S. (2000). On the computation of high order pseudospectral derivatives. Applied Numerical Mathematics, 33(1), 151-159. [ Links ]

Denlinger, R. P. & Iverson, R. M. (2001). Flow of variably fluidized granular masses across three-dimensional terrain: Numerical predictions and experimental tests. Journal of Geophysical Research, 106(B1), 553-566. [ Links ]

Denlinger, R. P. & Iverson, R. M. (2004). Granular avalanches across irregular three-dimensional terrain: 1. theory and computation. Journal of Geophysical Research, 109 (F1). [ Links ]

Diamessis, P. J., Lin, Y.-C., & Domaradzki, J. A. (2008). Effective numerical viscosity in spectral multidomain penalty method-based simulations of localized turbulence. Journal of Computational Physics, 227(17), 8145-8164. [ Links ]

Diamessis, P. J. & Redekopp, L. G. (2006). Numerical investigation of solitary internal wave-induced global instability in shallow water benthic boundary layers. Journal of physical oceanography, 36(5), 784-812. [ Links ]

Escobar-Vargas, J., Diamessis, P., & Giraldo, F. (2012). High-order discontinuous element-based schemes for the inviscid shallow water equations: Spectral multidomain penalty and discontinuous Galerkin methods. Applied Mathematics and Computation, 218(9), 4825-4848. [ Links ]

Fletcher, C. (1991). Computational techniques for fluid dynamics, volume 1 of Springer series in computational physics. Springer-Verlag, Berlin; New York, 2nd edition, 409 pp. [ Links ]

Giraldo, F. & Warburton, T. (2008). A high-order triangular discontinuous Galerkin oceanic shallow water model. International Journal for Numerical Methods in Fluids, 56(7), 899-925. [ Links ]

Gottlieb, D., & Hesthaven, J. S. (2001). Spectral methods for hyperbolic problems. Journal of Computational and Applied Mathematics, 128(1- 2), 83-131. [ Links ]

Gottlieb, D., & Shu, C. W. (1997). On the Gibbs phenomenon and its resolution. SIAM review, 39(4), 644-668. [ Links ]

Gottlieb, D., & Tadmor, E. (1985). Recovering pointwise values ofdiscontinuous data within spectral accuracy. In Progress and supercomputing in computational fluid dynamics, 357-375. [ Links ]

Gray, J., Wieland, M., & Hutter, K. (1999). Gravity-driven free surface flow of granular avalanches over complex basal topography. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 455(1985), 93-138. [ Links ]

Gray, J. M. N. T. (2003). Rapid granular avalanches. In Dynamic Response of Granular and Porous Materials under Large and Catastrophic Deformations, 3-42. [ Links ]

Springer. Hesthaven, J. S. & Gottlieb, A. D. (1996). A stable penalty method for the compressible navier-stokes equations: I. open boundary conditions. Journal on Scientific Computing, 17(3), 579-612. [ Links ]

Hesthaven, J. S., Gottlieb, S., & Gottlieb, D. (2007). Spectral methods for time-dependent problems (Vol. 21). Cambridge University Press, 272 pp. [ Links ]

Hesthaven, J. S. & Kirby, R. (2008). Filtering in Legendre spectral methods. Mathematics of Computation, 77(263), 1425-1452. [ Links ]

Hesthaven, J.S.& Warburton, T. (2007). Nodal discontinuous Galerkin methods: algorithms, analysis, and applications. Springer Science & Business Media. 495 pp. [ Links ]

Hutter, K. & Greve, R. (1993). Two-dimensional similarity solutions for finite-mass granular avalanches with coulomb-and viscous-type frictional resistance. Journal of Glaciology, 39(132), 357-372. [ Links ]

Hutter, K. & Koch, T. (1991). Motion of a granular avalanche in an exponentially curved chute: experiments and theoretical predictions. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 334(1633), 93-138. [ Links ]

Hutter, K.,Wang, Y., & Pudasaini, S. P. (2005). The Savage-Hutter avalanche model: how far can it be pushed? Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 363(1832), 1507-1528. [ Links ]

Iverson, R. M. (2014). Debris flows: behaviour and hazard assessment. Geology Today, 30(1), 15-20. [ Links ]

Iverson, R. M., & Denlinger, R. P. (2001). Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory. Journal of Geophysical Research: Solid Earth(1978-2012), 106(B1), 537-552. [ Links ]

Koch, T., Greve, R., & Hutter, K. (1994). Unconfined flow of granular avalanches along a partly curved surface. II. Experiments and numerical computations. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 445(1924), 415-435. [ Links ]

Koschdon, K. & Schafer, M. (2003). A Lagrangian-Eulerian finite-volume method for simulating free surface flows of granular avalanches. In Dynamic Response of Granular and Porous Materials under Large and Catastrophic Deformations, 83-108. Springer. [ Links ]

Kundu, P. & Cohen, L. (2008). Fluid Mechanics. Elsevier Acad. Press, London, England, 729 pp. [ Links ]

LeVeque, R. J. (2002). Finite volume methods for hyperbolic problems, volumen 31. Cambridge university press, Cambridge, England, 529 pp. [ Links ]

Maday, Y., Kaber, S. M. O., & Tadmor, E. (1993). Legendre pseudo spectral viscosity method for nonlinear conservation laws. SIAM Journal on Numerical Analysis, 30(2), 321-342. [ Links ]

Navarra, A., Stern, W., & Miyakoda, K. (1994). Reduction of the Gibbs oscillation in spectral model simulations. Journal of climate, 7(8), 1169-1183. [ Links ]

Ouyang, C., He, S., Xu, Q., Luo, Y., & Zhang, W. (2013). A MacCormack-TVD finite difference method to simulate the mass ow in mountainous terrain with variable computational domain. Computers & Geosciences, 52(0), 1-10. [ Links ]

Pastor, M., Blanc, T., Haddad, B., Drempetic, V., Morles, M. S., Dutto, P., Stickle, M. M., Mira, P., & Merodo, J. F. (2015). Depth averaged models for fast landslide propagation: mathematical, rheological and numerical aspects. Archives of Computational Methods in Engineering, 22(1), 67-104. [ Links ]

Pudasaini, S. P. (2012). A general two-phase debris ow model. Journal of Geophysical Research: Earth Surface, 117(F3), 799-819. [ Links ]

Pudasaini, S. P. & Hutter, K. (2007). Avalanche dynamics: dynamics of rapid flows of dense granular avalanches. Springer, Berlin, 571 pp. [ Links ]

Reed, W. H. & Hill, T. (1973). Triangular mesh methods for the neutron transport equation. Technical report, Los Alamos Scientific Lab., N. Mex. (USA) . [ Links ]

Savage, S. B. & Hutter, K. (1989). The motion of a finite mass of granular material down a rough incline. Journal of Fluid Mechanics, 199, 177-215. [ Links ]

Spiteri, R. J. & Ruuth, S. J. (2002). A new class of optimal high order strong stability-preserving time discretization methods. SIAM Journal on Numerical Analysis, 40(0), 469-491. [ Links ]

Tai, Y. & Gray, J. (1998). Limiting stress states in granular avalanches. Annals of Glaciology, (26), 272-276. [ Links ]

Tai, Y.-C., Hutter, K., & Gray, J. (2001). Dense granular avalanches: mathematical description and experimental validation. In Geomorphological fluid mechanics, 339-366. Springer. [ Links ]

Tai, Y.-C., Noelle, S., Gray, J., & Hutter, K. (2002). Shock-capturing and front-tracking methods for granular avalanches. Journal of Computational Physics, 175(1), 269-301. [ Links ]

Tai, Y.-C., Wang, Y., Gray, J., & Hutter, K. (1999). Methods of similitude in granular avalanche flows. In Advances in Cold-Region Thermal Engineering and Sciences, 415-428. [ Links ]

Toro, E. F. (2013). Riemann solvers and numerical methods for fluid dynamics: a practical introduction. Springer Science & Business Media, Dordrecht; New York, 3rd edition. [ Links ]

Trujillo-Vela, M. G. (2015). Simulación numérica de flujos de material desagregado. Método penalizado de multidominios espectrales. M.Sc. Thesis, Maestría en Hidrosistemas, Pontificia Universidad Javeriana, Bogotá D. C., Colombia [ Links ]

Vandeven, H. (1991). Family of spectral filters for discontinuous problems. Journal of Scientific Computing, 6(2), 159-192. [ Links ]

Vollmoller, P. (2004). A shock-capturing wave-propagation method for dry and saturated granular flows. Journal of Computational Physics, 199(1), 150-174. [ Links ]

Wang, Y., Hutter, K., & Pudasaini, S. P. (2004). The Savage-Hutter theory: A system of partial differential equations for avalanche flows of snow, debris, and mud. Zeitschrift fur angewandte Mathematik und Mechanik ZAMM, 84(8), 507-527. [ Links ]

Wieland, M., Gray, J., & Hutter, K. (1999). Channelized free-surface flow of cohesionless granular avalanches in a chute with shallow lateral curvature. Journal of"Fluid Mechanics, 392, 73-100. [ Links ]

Zhou, J. G. (2004). Lattice Boltzmann methods for shallow water flows. Springer, New York. [ Links ]

How to cite item Trujillo-Vela, M. G., Escobar-Vargas, J. A., & Ramos-Cañón, A. M. (2019). A spectral multidomain penalty method solver for the numerical simulation of granular avalanches. Earth Sciences Research Journal, 23(4), 317329. DOI: https://doi.org/10.15446/esrj.v23n4.77683

Appendix A. Spectral multidomain penalty method

In order to assess the potential of high-order discontinuous methods for real applications and in order to reduce the numerical error introduced by low-order techniques in this paper, we performed the solution of the granular flow equations (equation 1) with a method that combines an exponential convergence, a weak artificial dissipation, and an adaptability based in classic finite elements-volumes techniques (Boyd, 2001). The spectral multidomain penalty method (SMPM) 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 (Hesthaven and Gottlieb, 1996, Escobar-Vargas et al., 2012). The formulation of SMPM will be presented in two parts: the subdomain interior and the interfacial treatments.

Subdomain interior. In the interior of each subdomain any function q(x, y, t) can be approximated by using N-th order Lagrange interpolating polynomials on a Gauss-Lobatto-Legendre grid (GLL) (Hesthaven and Gottlieb, 1996).

where q(x i, y j, t) is the value of function at the discrete point x i, y j ), and li(x), l j (y) are the i-th and j-th Lagrange interpolating polynomials based on the GLL nodes in the and y directions, respectively. The spatial derivatives in the x-direction in the global coordinate system are approximated as

where ξ∂/∂x represents the mapping from the local coordinates system I ξЄ [-1,1], given by GLL points into the global coordinates system x Є ℝ.D ik is the Legendre spectral differentiation matrix, which is computed following the procedure suggested by Costa and Don (2000).

Interfacial treatment. The penalized form of the granular flow equations in a collocation point situated on a subdomain boundary requires an additional penalization term, which allows us establishes the continuity between adjacent points (Escobar-Vargas, 2012).

where n (e,l ) is the unit vector at the subdomain boundary, l represents each one of the four edge of the domain, represents the information flux between adjacent points, and is the penalization coefficient (Escobar- Vargas, 2012).

where ω are weights associated to the GLL points. For a uniform polynomial approximation, N, in each subdomain can be used a unique value of ω = 2/N{N +1) and ∆s = {∆x, ∆y) which depends on the orientation of the interfaces in the subdomains.

Appendix B. Strong stability preserving Runge-Kutta

In order to maintain the precision offered by the spatial discretization method, it is used a strong stability preserving Runge-Kutta method (SSP-RK) for the temporal discretization. By considering the following initial value problem

In this scheme, the prediction at the time t + 1 (i.e., q t+1) is based on both the existing solution at the time t (i.e., qt) and the evaluation of the forcing terms R(q) at different stages (Spiteri and Ruuth, 2002).

For consistency, we must have that = 1, i= 1, 2,..., s, where s is the number of stages of the SSP-RK approach, α¡,k and β ¡,k are constant coefficients where all α¡,k ≥ 0 and α¡,k = 0 only if β, ¡,k = 0. Spiteri and Ruuth (2002) give these coefficients. ∆t is the size of the time-step.

Appendix C. Stabilization technique

There are multiples methods for improving accuracy and stability of spectral and pseudo spectral techniques, such as physical space filters (Gottlieb and Tadmor, 1985, Gottlieb and Shu, 1997), spectral vanishing viscosity methods (Maday et al., 1993), dealiasing (Canuto et al., 2012) and over integration (Hesthaven and Kirby, 2008). A family of regularization techniques related to the first two above-mentioned is the spectral filters. Since there is not a dissipative term in the SMPM to remove the artificial oscillations, a spectral filter used by Vandeven, (1991), Gottlieb and Hesthaven, (2001), Diamessis et al., (2008) is necessary to reduce the high-frequency modes that appear close to the discontinuities (Hesthaven and Warburton, 2007). It was selected this spectral filter (exponential) because of the stability and robustness that it adds to the method (Hesthaven and Kirby, 2008, Hesthaven et al., 2007). A comparison of different spectral filtering methods and discussion of the advantages and weaknesses of using them is performed in Hesthaven and Kirby, (2008), Blackburn and Schmidt, (2003), Navarra et al., (1994).

Based on the results of this work, it can be concluded that the spectral filter allows us to obtain stable solutions, whereas few solutions are obtained when a filter is not used. Thus, it is possible to avoid artificial and persistent oscillations around points where discontinuities appear. Also, the error in mass conservation has a smoother behavior over time when using a filter.

Received: February 05, 2019; Accepted: September 20, 2019

*Corresponding author: mario.tmjillo@javeriana.edu.co

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License