Introduction

The Savage-Hutter equations (1989, 1991) are a set of hyperbolic evolution equations that describe the geometry and depth-averaged velocity in finite-mass avalanche flows. The original equations are not written in conservative form and thus can also only reproduce smooth solutions in cases when shock should form. All numerical schemes based on these non-conservative forms generated high oscillations at numerical instabilities in the vicinity where shocks were expected. These were at first damped by adding numerical diffusion but not basing the model on equations in conservative form. This has now been done, and the numerical integration scheme allows the shocks to be identified without the manifestation of spurious oscillations or the formation of instabilities and without the introduction of additional numerical diffusion. As a result, the numerically determined solutions reflect the properties of the model without diffusive contamination.

Shocks are initiated when the avalanche velocity is faster than its characteristic speed and the avalanche front reaches the base of the slope or a solid wall. The shock wave then travels upslope through the avalanche. Many detailed investigations of granular shocks were made by Reference Gray and HutterGray and Hutter (1997), in which the shock waves are considered to be an important feature in the granular flows. The Savage-Hutter equations of granular avalanche flows are a hyperbolic system of equations, such that velocities may be supercritical and shock waves may occur.

To avoid the numerical instability, an artificial viscosity term μ2u/x2 is introduced in the Lagrangian finite-difference method proposed by Reference Savage and HutterSavage and Hutter (1989, Reference Savage and Hutter1991) for numerical stability, where the artificial viscosity μ was found to have values of 0.01–0.03. Such instability is here supposed to be caused by shocks. Guido and Bartelt (unpublished) applied an upwind total variation diminishing (TVD) difference scheme to solve quasi-Savage-Hutter equations and concluded the numerical stability without using artificial damping with their discretization schemes with respect to the characteristic variables. However, in the quasi-Savage-Hutter and Savage-Hutter equations there are jumps within the characteristics as well as in the characteristic variables. The treatment was not clear in their statement. An alternative is to solve the equations with the conservative variables, the thickness h and the thickness-integrated moment m = hu. In this paper a non-oscillatory central (NOC) shock -capturing numerical scheme (Reference Nessyahu and TadmorNessyahu and Tadmor, 1990) is applied to the Savage-Hutter equations (1989, 1991) with respect to the conservative variables to describe the shock formation directly without adding artificial viscosity or other regularizations. Since this numerical scheme is based on a stationary uniform mesh, it is impossible to point out the margin location (the material-free region) without extra treatment. Therefore we combine our NOC scheme with a front-tracking (FT) method, also developed in Reference TaiTai (2000) and Reference Tai, Noelle, Gray and HutterTai and others (2000). The resulting NOC-FT scheme is able to accurately determine the material-free region.

Governing Equations

In the Savage-Hutter theory (1989, 1991) the dry cohesionless granular material is assumed to be an incompressible continuum with constant density throughout the entire body. During flow the body behaves as a Mohr-Coulomb plastic material at yield and slides over a rigid basal topography with similar Coulomb-type frictional behaviour. Scaling analysis isolates the physically significant terms in the governing equations and identifies those terms that can be neglected. Finally, depth integration reduces the theory to one spatial dimension.

The leading-order, dimensionless, depth-integrated mass balance reduces to

where h is the avalanche thickness and u is the velocity in the downslope direction. The leading-order, dimensionless, depth-integrated momentum balance is

with net driving force

where ζ, and λK are local slope inclination angle, basal Coulomb dry-friction angle and local curvature of the chute profile, respectively. The term sgn(u) comes from the dry Coulomb drag friction, and ε≪ 1 is the aspect ratio of the typical thickness and the length of the avalanche. Note that Equations (1) and (2) are written in conservative form, while in the original Savage-Hutter theory the smoothness assumption yields the momentum-balance equation to the velocity-evolution equation

The factor β_x is defined as βx = ε cos ζKx, and the earth pressure coefficient K_x relates the in-plane stresses to the normal stress within the avalanche K_x = p_xx/p_zz. The value of K_x is given by the Mohr-Coulomb criterion with an ad hoc assumption, depending on whether the avalanche body extends or contracts,

with

where φ is the internal friction angle of the granular material.

The original Lagrange finite-difference scheme (Reference Savage and HutterSavage and Hutter, 1989) is implemented by the equation system, (1) and (4), in Lagrangian form with primitive variables, i.e. the local thickness of the avalanche h, and the downslope velocity u. However, the shock-capturing scheme employed here is applied to this equation system in conservative form, (Equations (1) and (2)), with conservative variables, i.e. the thickness h and the thickness-integrated moment m = hu.

Shock-Capturing Numerical Scheme

To describe the shock formation the NOC differencing scheme (Reference Nessyahu and TadmorNessyahu and Tadmor, 1990) is here employed, which is similar in form to the Lax-Friedrichs scheme. For completeness, we briefly review the scheme. The fundamental difference between the NOC scheme and standard upwind finite-difference schemes is that the boundaries of the cells at the new time level are the centers of the cells at the old time level (see Fig. 1). In other words, the cell average values bounded by [xj, xj+i at time level tn+1 are computed from the information from the cell average values bounded by bounded by [xj+l/2, xj+3/2] at time level tn (see Fig. 1, top). This motivates the term “central”. For the NOC scheme, the expensive Riemann solvers used in upwind schemes can be completely avoided. The resulting scheme is easy to code, computationally efficient and can be applied to general systems of conservation laws, where the solution of the Riemann problem (i.e. the initial value problem with piecewise constant data) may be complicated or even impossible. For further details and references to recent related work we refer the reader to Reference Tai, Noelle, Gray and HutterTai and others (2000).

Fig. 1.

Diagram of NOC scheme. ToP: Gridpoints computed in the NOC method. Bottom: NOC scheme computational diagram, where • indicate the gridpoints at time level n and n H represent the positions where the fluxes fat time level n +1/2 are approximated, ◊ are for the source term s, and A denote the quarter and three-quarter points meaning (e.g.

For simplicity the Savage-Hutter equations in the form of Equations (1) and (2) are considered with w = (h, m)_T as basic variables and written in the form

where

Let denote the cell average oft.; at (xj, t1); integrating the governing equations over the rectangle [x3,xJ+1]× [tn,tn+1] gives

So, the cell average over can be computed by the following discretized equation:

as illustrated in Figure 1 (bottom). The values of and are determined by the reconstruction over the jth and (j + 1)th cell,i.e.

where denotes the cell mean derivative determined by TVD hmiter (e.g. Reference YeeYee, 1989) or the weighted essentially non-oscillatory (WENO) cell reconstruction (Reference Jiang and ShuJiang and Shu, 1996).

The transport flux f is approximated by the physical quantities at (xj, tn+1/2) i.e.

The integration of the source term represents the volume of s over [xj, xj=1] which is approximated by the values at

and

The temporal derivative in Equations (12) and (14) is determined by using Equation (7),

where

To keep the numerical stability and to ensure the non-oscillatory property at a discontinuity the following CFL condition

must hold during computation, where

is the maximum wave speed.

Front-Tracking

The shock-capturing method discretizes the governing equations on a stationary uniform mesh. The margins may be located between gridpoints, so that it is impossible to point out the margin location without extra treatment. In the front-tracking method used in this paper and developed inReference Tai, Noelle, Gray and HutterTai and others (2000), the margin cell is reconstructed by a piecewise linear distribution, so that outside the avalanche body there is no material and the cell average does not change due to the reconstruction (e.g. see Fig. 2 for the front margin). It can be easily proven that the margin point fulfils the smoothness assumption (Reference TaiTai, 2000), and consequently its motion can be described by the following evolution equation

which is equivalent to the momentum-balance Equation (4) in the original Savage and Hutter theory. This combination is then called the NOG-FT method. During computation all the numerical fluxes and the source term must be considered if the margin moves into the neighbouring cell.

Fig. 2.

Cell reconstrution of the depth hf(x) within the front-margin fth cell. In the front-tracking method a piecewise linear reconstrution hAx) over the margin cell is well defined, such that the value of the cell average does not change due to the definition of the reconstruction.

Numerical Simulations

Conclusions

The shock-capturing numerical scheme provides an excellent tool to describe shock formation in granular flows with the Savage and Hutter theory. Combining this numerical scheme with the front-tracking method makes it possible to determine the margin locations of the avalanche. Numerical tests of the upward-moving shock wave and a comparison with the parabolic similarity solution confirm the superlority of this combination over previous schemes. Further numerical experiments, including flows for which the basal and internal friction angles are different from each other, are reported inReference TaiTai (2000) and Reference Tai, Noelle, Gray and HutterTai and others (2000).