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 Kx relates the in-plane stresses to the normal stress within the avalanche Kx = pxx/pzz. The value of Kx 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
Upward-moving shock wave
Shock formations are often observed when the avalanche slides into the run-out horizontal zone. Here the front part comes to rest, while the tail accelerates further and its velocity becomes supercritical. A comparison is made between the shock-capturing method and the Lagrangian moving grid technique. While such shock formation normally induces numerical instability for the Lagrangian method, the shock-capturing scheme can behave well for such shocks.
The granular material released from a parabolic cap slides down an inclined plane and merges into the run-out horizontal zone. The centre of the cap is initially located at x = 4.0, and the initial radius and the height are 3.2 and 2.2 dimensionless length units, respectively. The inclination angle of the inclined plane is 40°, and the transition region lies between x = 21.5 and x = 25.5. The basal and internal friction angle are both 30°.
Figure 3 illustrates the simulated process as the avalanche slides on the inclined plane into the horizontal runout zone. The avalanche body extends on the inclined plane until the front reaches the horizontal run-out zone. Here the basal friction is enough to bring the front of the granular material to rest, but the part of the rear accelerates further. At this stage, a shock (surge) wave is created (t = 12), which moves upward. Such shock waves make the Lagrangian method unstable if no artificial viscosity is applied (see the lowermost panel in Fig. 3).
Fig. 3. Process of the avalanche simulated by the shock-capturing NOC method at t = 0, 4, 8, 12, 16, 20 dimensionless time units. As the front reaches the run-out zone and comes to rest, the part of the tail accelerates further and the avalanche body contracts. Once the velocity becomes supercritical, a shock wave develops, which moves upward. For the Lagrangian moving-grid method without applying any artificial viscosity, oscillations take place (lowermost panel) when the shock occurs, which may result in an element boundary overtaking the other one and destroying the scheme.
Parabolic similarity solution
To test our shock-capturing method for a smooth solution with moving boundaries, the results are compared with the parabolic similarity solution (Reference Savage and HutterSavage and Hutter, 1989), which has already been successfully modelled using the Lagrangian technique (Reference Savage and HutterSavage and Hutter, 1989,Reference Savage and Hutter1991). We refer to Reference TaiTai (2000) for a generalization of that similarity solution.
Consider a finite mass of granular avalanche sliding on an infinitely long inclined plane at C = 40°, where the internal and basal angles of friction are φ = δ = 30°. With a linearly increasing distribution of velocity, ∂u/∂x > 0, the avalanche body preserves its parabolic form.
Figure 4 demonstrates the results computed by the NOC scheme without the front-tracking method, where the TVD superbee (e.g. Reference YeeYee, 1989) limiter is employed and a thin layer ho = 10–4 is added to the entire computational domain. The results show that the added thin layer does not influence the depth profile very much, if it is sufficiently small, but one cannot determine the margin locations.
Fig. 4. Depth (upper) and velocity (lower) profiles of the parabolic similarity solution computed by the NOC scheme with the TVD Superbee limiter, where the whole computational domain is divided into 90 cells, the circles denote simulated results and the solid lines represent the exact solution. A thin layer with h0 = 10–4 is added to the whole computational domain.
However, as can be seen in the velocity profiles of Figure 4, there are large velocity gradients around the margins. Since the regions outside the margin are covered by the uniformly added thin layer, there is no contribution from the depth gradient, ∂h/∂x, in the momentum-balance Equation (2). On the other hand, inside the avalanche body there is a permanent contribution from the depth gradient for the momentum balance equation. The velocity in this region will therefore tend to increase differently from that outside the avalanche body, so a jump of velocity develops around the margins. Here the jumps are reproduced by a smeared transition, (see the velocity profiles in Fig. 4). This is why there are deviations of the depth around the margin. Furthermore, the smeared transition of velocity results in several cells with ∂u/∂x < 0 around the margins. This violates the presumption ∂u/∂x > 0 in the parabolic similarity solution problem. Therefore, the front-tracking method must be introduced.
Near the centre of the avalanche body there is a slight buckling. This is caused by the first-order accurate approximation at the local extremum and is called the clipping phenomenon, which is the weakpoint of the TVD schemes. One way of avoiding this weakness is to apply the higher-order WENO cell reconstruction.
Figure 5 demonstrates the results (circles) from the NOG scheme with the front-tracking method, where a WENO piecewise quadratic interpolation (Reference Jiang and ShuJiang and Shu, 1996) is employed. With this method both the parabolic form of the avalanche depth and the margin locations can be well described. The so-called clipping phenomenon is successfully removed. The presumed stress state, ∂u/∂x > 0, in the parabolic similarity solution problem is also maintained.
Fig. 5. Depth (upper) and velocity (lower) profiles of the parabolic similarity solution at t = 6 computed by the NOC scheme with front-tracking method, where the piecewise quadratic (WENO r = 3) interpolations are implemented. The circles denote simulated results, and the solid lines represent the exact solution. The whole computational domain is divided into 90 cells, and the Courant number is selected to be 0.3.
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).
Acknowledgements
This research was supported by the Deutsche Forschungsgemeinschaft projects SFB 298 “Deformation und Versagen bei metalhschen und granularen Strukturen” at Darmstadt University of Technology and SFB 256 “Nichthneare Partielle Differentialgleichungen” at Bonn University
