A Finite Volume Scheme for Savage-Hutter Equations on Unstructured Grids

TL;DR

This paper introduces a finite volume numerical scheme on unstructured grids for Savage-Hutter equations, effectively handling discontinuities and non-Galilean invariance, demonstrated through granular avalanche flow simulations.

Contribution

A novel Godunov-type finite volume scheme with a modified HLL flux for Savage-Hutter equations on unstructured grids, addressing non-Galilean invariance and discontinuities.

Findings

Scheme accurately captures shock waves in granular flows

Numerical results agree well with reference solutions

Applicable to complex avalanche simulations

Abstract

A Godunov-type finite volume scheme on unstructured triangular grids is proposed to numerically solve the Savage-Hutter equations in curvilinear coordinate. We show the direct observation that the model is a not Galilean invariant system. At the cell boundary, the modified Harten-Lax-van Leer (HLL) approximate Riemann solver is adopted to calculate the numerical flux. The modified HLL flux is not troubled by the lack of Galilean invariance of the model and it is helpful to handle discontinuities at free interface. Rigidly the system is not always a hyperbolic system due to the dependence of flux on the velocity gradient. Even though, our numerical results still show quite good agreements to reference solutions. The simulations for granular avalanche flows with shock waves indicate that the scheme is applicable.

Figures

Equations

$$\frac{\partial h}{\partial t}+\frac{\partial }{\partial x}(hu)+\frac{\partial }{\partial y}(hv)=0,$$

$$\frac{\partial (hu)}{\partial t}+\frac{\partial }{\partial x}(h{u}^2+\frac{1}{2}{\beta }_x{h}^2)+\frac{\partial }{\partial y}(huv)=h{s}_x,$$

$$\frac{\partial (hv)}{\partial t}+\frac{\partial }{\partial x}(huv)+\frac{\partial }{\partial y}(h{v}^2+\frac{1}{2}{\beta }_y{h}^2)=h{s}_y,$$

$${\beta }_x=ϵ\cos \zeta K_x,{\beta }_y=ϵ\cos \zeta K_y,$$

$$K_{x_{act/pass}}=2\left(1\mp \sqrt{1-\frac{{\cos }^2\varphi }{{\cos }^2\delta }}\right){\sec }^2\varphi -1,$$

$$K_{y_{act/pass}}=\frac{1}{2}\left(K_x+1\mp \sqrt{(K_x-1{)}^2+4{\tan }^2\delta }\right).$$

$$K_{x}=\left\{\begin{array}[]{lr}K_{x_{act}},&\partial u/\partial x\geq 0,\\ K_{x_{pass}},&\partial u/\partial x<0,\end{array}\right.$$

$$K_{y}=\left\{\begin{array}[]{llr}K_{y_{act}}^{x_{act}},&\partial u/\partial x\geq 0,&\partial v/\partial y\geq 0,\\[5.69054pt] K_{y_{pass}}^{x_{act}},&\partial u/\partial x\geq 0,&\partial v/\partial y<0,\\[5.69054pt] K_{y_{act}}^{x_{pass}},&\partial u/\partial x<0,&\partial v/\partial y\geq 0,\\[5.69054pt] K_{y_{pass}}^{x_{pass}},&\partial u/\partial x<0,&\partial v/\partial y<0.\end{array}\right.$$

$$s_x=\sin \zeta -\frac{u}{|\mathbf{u}|}\tan \delta (\cos \zeta +\lambda \kappa u^2)-ϵ\cos \zeta \frac{\partial z_b}{\partial x},$$

$$s_y=-\frac{v}{|\mathbf{u}|}\tan \delta (\cos \zeta +\lambda \kappa u^2)-ϵ\cos \zeta \frac{\partial z_b}{\partial y},$$

$$\frac{\partial \mathbf{U}}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{U})}{\partial x}+\frac{\partial \mathbf{G}(\mathbf{U})}{\partial y}=\mathbf{S}(\mathbf{U}),$$

$$\displaystyle\bm{U}=\left[\begin{array}[]{c}h\\ hu\\ hv\end{array}\right],\quad\quad\bm{F}=\left[\begin{array}[]{c}hu\\ hu^{2}+\frac{1}{2}\beta_{x}h^{2}\\ huv\end{array}\right],$$

$$\displaystyle\bm{G}=\left[\begin{array}[]{c}hv\\ huv\\ hv^{2}+\frac{1}{2}\beta_{y}h^{2}\end{array}\right],\quad\quad\bm{S}=\left[\begin{array}[]{c}0\\ hs_{x}\\ hs_{y}\end{array}\right].$$

$$\bm{J}=\frac{\partial(\bm{F}n_{x}+\bm{G}n_{y})}{\partial\bm{U}}=\left[\begin{array}[]{ccc}0&n_{x}&n_{y}\\ (\beta_{x}h-u^{2})n_{x}-uvn_{y}&2un_{x}+vn_{y}&un_{y}\\ -uvn_{x}+(\beta_{y}h-v^{2})n_{y}&vn_{x}&un_{x}+2vn_{y}\end{array}\right],$$

$${\lambda }_1$$

$${\lambda }_2$$

$${\lambda }_3$$

$$\cos \theta \mathbf{F}(\mathbf{U})+\sin \theta \mathbf{G}(\mathbf{U})={\mathbf{T}}^{-1}\mathbf{F}(\mathbf{TU})$$

$$\bm{T}=\left[\begin{array}[]{ccc}1&0&0\\ 0&\cos\theta&\sin\theta\\ 0&-\sin\theta&\cos\theta\end{array}\right].$$

$$\bm{TU}=\left[\begin{array}[]{c}h\\ h(u\cos{\theta}+v\sin{\theta})\\ h(v\cos{\theta}-u\sin{\theta})\end{array}\right].$$

$$\bm{F}(\bm{TU})=\left[\begin{array}[]{c}h(u\cos{\theta}+v\sin{\theta})\\ \frac{h^{2}\beta_{x}}{2}+\frac{(hu\cos{\theta}+hv\sin{\theta})^{2}}{h}\\ \frac{(hv\cos{\theta}-hu\sin{\theta})(hu\cos{\theta}+hv\sin{\theta})}{h}\end{array}\right].$$

$$\bm{T^{-1}F}(\bm{TU})=\left[\begin{array}[]{c}h(u\cos{\theta}+v\sin{\theta})\\ \frac{1}{2}h(2u^{2}+\beta_{x})\cos{\theta}+huv\sin{\theta}\\ huv\cos{\theta}+\frac{1}{2}h(2v^{2}+h\beta_{x})\sin{\theta}\end{array}\right].$$

$$\cos\theta\bm{F}(\bm{U})+\sin\theta\bm{G}(\bm{U})=\left[\begin{array}[]{c}h(u\cos{\theta}+v\sin{\theta})\\ \frac{1}{2}h(2u^{2}+\beta_{x})\cos{\theta}+huv\sin{\theta}\\ huv\cos{\theta}+\frac{1}{2}h(2v^{2}+h\beta_{y})\sin{\theta}\end{array}\right].$$

$${{\mathbf{U}}^n|}_{{\tau }_i}(\mathbf{x})={\mathbf{U}}_i^n+\nabla {\mathbf{U}}_i^n\left(\mathbf{x}-{\mathbf{x}}_i\right)$$

$$\frac{\partial }{\partial t}{\int }_{{\tau }_i}\mathbf{U}d\Omega +{\int }_{{\tau }_i}\left(\frac{\partial \mathbf{F}}{\partial x}+\frac{\partial \mathbf{G}}{\partial y}\right)d\Omega ={\int }_{{\tau }_i}\mathbf{S}d\Omega .$$

$$\frac{\partial }{\partial t}{\int }_{{\tau }_i}\mathbf{U}d\Omega +{\int }_{\partial {\tau }_i}\mathcal{F}\cdot \mathbf{n}ds={\int }_{{\tau }_i}\mathbf{S}d\Omega ,$$

$${\mathbf{U}}_i=\frac{1}{|{\tau }_i|}{\int }_{{\tau }_i}\mathbf{U}d\Omega ,$$

$${\mathbf{U}}_i^{n+1/2}={\mathbf{U}}_i-\frac{\Delta t}{2|{\tau }_i|}\sum _{j=1}^3{\int }_{\partial {\tau }_i,j}\mathcal{F}({\mathbf{U}}_{in}^n)\cdot {\mathbf{n}}_{j}ds+\frac{\Delta t}{2}{\mathbf{S}}_i^n,$$

$${\mathbf{U}}_i^{n+1}={\mathbf{U}}_i^n-\frac{\Delta t}{|{\tau }_i|}\sum _{j=1}^3{\int }_{\partial {\tau }_i,j}{\mathcal{F}}_j^{\ast }({\mathbf{U}}_{in}^{n+1/2},{\mathbf{U}}_{out}^{n+1/2})\cdot {\mathbf{n}}_{j}ds+\Delta t{\mathbf{S}}_i^{n+1/2}.$$

$$\bm{\mathcal{F}}^{*}\cdot\bm{n}=\left\{\begin{array}[]{lc}\bm{\mathcal{F}}(\bm{U}_{L})\cdot\bm{n},&\mathrm{if}\quad s_{L}\geq 0,\\ \dfrac{s_{R}\bm{\mathcal{F}}(\bm{U}_{L})\cdot\bm{n}-s_{L}\bm{\mathcal{F}}(\bm{U}_{R})\cdot\bm{n}+s_{R}s_{L}(\bm{U}_{R}-\bm{U}_{L})}{s_{R}-s_{L}},&\mathrm{if}\quad s_{L}\leq 0\leq s_{R},\\ \bm{\mathcal{F}}(\bm{U}_{R})\cdot\bm{n},&\mathrm{if}\quad s_{R}\leq 0,\end{array}\right.$$

$$s_{L}=\left\{\begin{array}[]{lc}\min(\bm{u}_{L}\cdot\bm{n}-\sqrt{c_{L}h_{L}},u_{*}-\sqrt{c_{L}h_{*}}),&\mathrm{if}\quad h_{L},h_{R}>0,\\ \bm{u}_{L}\cdot\bm{n}-\sqrt{c_{L}h_{L}},&\mathrm{if}\quad h_{R}=0,\\ \bm{u}_{R}\cdot\bm{n}-2\sqrt{c_{R}h_{R}},&\mathrm{if}\quad h_{L}=0,\end{array}\right.$$

$$s_{R}=\left\{\begin{array}[]{lc}\max(\bm{u}_{R}\cdot\bm{n}+\sqrt{c_{R}h_{R}},u_{*}+\sqrt{c_{R}h_{*}}),&\mathrm{if}\quad h_{L},h_{R}>0,\\[5.69054pt] \bm{u}_{L}\cdot\bm{n}+2\sqrt{c_{L}h_{L}},&\mathrm{if}\quad h_{R}=0,\\[5.69054pt] \bm{u}_{R}\cdot\bm{n}+\sqrt{c_{R}h_{R}},&\mathrm{if}\quad h_{L}=0,\end{array}\right.$$

Full text

keywords: granular avalanche flow, Savage-Hutter equations, finite volume method, Galilean invariant

1. Introduction

The Savage-Hutter model was first proposed to describe the motion of a finite mass of granular material flowing down a rough incline, in which the granular material was treated as an incompressible continuum [25, 7, 11]. It was then extended to 2D and with curvilinear coordinate [10, 9, 8] which may be applied to study natural disasters as landslides and debris flows [18]. For more details on Savage-Hutter model, we refer to [11, 23, 6, 12, 13, 22].

Numerical method for the Savage-Hutter model may be dated back to work in [25, 15] using finite difference methods which are not able to capture shock waves. Later on, regardless of the Savage-Hutter model or other avalanches of granular flow models, the Godunov-type schemes are commonly used in literature [4, 29, 1, 2, 30, 28, 24] considering that most of them have hyperbolic nature of the depth averaged equations. We note that rigidly, the Savage-Hutter model is not hyperbolic since its flux is depended on the gradient of velocity. Precisely, the dependence is on the sign of the gradient of velocity components, that the model does be hyperbolic in the region with monotonic velocity components. Some related interesting development on numerical methods can be found in such as [27, 19, 1, 20, 17].

In practical scenario, landslides, rock avalanches and debris flows usually occur in mountainous areas with complex topography, that we have to use unstructured grids for numerical simulation. However, there is few work of numerical methods for Savage-Hutter equations on unstructured grids. This may due primarily to the lack of Galilean invariance of the model, which is a direct observation to the model as we point out in section 2. In this paper, we are motivated to develop a high resolution and robust numerical model for Savage-Hutter equations under the unstructured grids. Whether the numerical method may provide a correct solution on the unstructured grids for a problem without Galilean invariance is a major question here. Meanwhile, the loss of hyperbolicity at the zeros of velocity gradients may some unpredictable behavior in numerical solutions, too.

Additionally, there are still some challenges in numerical solving strongly convective Savage-Hutter equations on unstructured grids. For example, shock formation is an essential mechanism in granular flows on an inclined surface merging into a horizontal run-out zone or encountering an obstacle when the velocity becomes subcritical from its supercritical state [29]. Therefore, numerical efficiency is an important issue to help resolve the steep gradients and moving fronts. Noticing that the Savage-Hutter equations and Shallow water equations have some similarities, we basically follow the method in [3] for shallow water equations. The modified HLL flux is adopted to calculate the numerical flux on the cell boundary. The formation of the modified HLL flux does not require the flux to be Galilean invariant, thus formally we have no difficulties in calculation of numerical flux. The techniques involved to improve efficiency include the MUSCL-Hancock scheme in time discretization, the ENO-type reconstruction and the $h$-adaptive method in spatial discretization.

We apply the numerical scheme to three examples. The first one is a granular dam break problem, and the second one is a granular avalanche flows down an incline plane and merges continuously into a horizontal plane. In the third example, a granular slides down an inclined plane merging into a horizontal run-out zone is simulated, whereby the flow is diverted by an obstacle which is located on the inclined plane. It is interesting that in these examples, the numerical solutions are agreed with the references solutions quite well. We are indicated that the scheme is applicable to the model, though the model is lack of the Galilean invariance and is not always hyperbolic. The reason why the loss of the hyperbolicity of the model are not destructive to the numerical methods is still under investigation.

The rest part of the paper is organized as follows. In section 2, the Savage-Hutter equations on the curivilinear coordinate is briefly introduced. We directly show that the system is not Galilean invariant. In section 3, we give the details of our numerical scheme. In section 4, we present the numerical results and a short conclusion in section 5 close the paper.

2. Governing Equations

There are various forms of Savage-Hutter equations according to different coordinate system and a detailed derivation of these different versions of Savage-Hutter equations has been given in [23]. Here, we confine ourselves to one of them. Precisely, the governing equations are built on an orthogonal curvilinear coordinate system (see Fig. 1) given in [1]. The curvilinear coordinate $Oxyz$ is defined on the reference surface, where the $x$-axis is oriented in the downslope, the $y$-axis lies in the cross-slope direction of the reference surface and the $z$-axis is normal to them. The downslope inclination angle of the reference surface $\zeta$ only depends on the downslope coordinates $x$, that is, there is no lateral variation in the $y$-direction, thus it is strictly horizontal. The shallow basal topography is defined by its elevation $z=z^{b}(x,y)$ above the curvilinear reference surface.

The corresponding non-dimensionless Savage-Hutter equations are

[TABLE]

[TABLE]

[TABLE]

where $h$ is the avalanche depth in the $z$ direction, and $\bm{u}=(u,v)$ are the depth-averaged velocity components in the downslope($x$) and cross-slope($y$) directions, respectively. The factors $\beta_{x}$, $\beta_{y}$ are defined as

[TABLE]

respectively, where $\epsilon$ is aspect ratio of the characteristic thickness to the characteristic downslope extent. Here $K_{x}$, $K_{y}$ are the $x$, $y$ directions earth pressure coefficients defined by the Mohr-Coulomb yield criterion. Generally, the earth pressure coefficients are considered in the active or passive state, which depends on whether the downslope and cross-slope flows are expanding or contracting [21]. Hutter et al. assumed that $K_{x}$, $K_{y}$ link the normal pressures in the $x$, $y$ directions with the overburden pressure and suggested that [26, 10]

[TABLE]

[TABLE]

where $\phi$ and $\delta$ are the internal and basal Coulomb friction angles, respectively. The subscripts “act” and “pass” denote active(“$-$”) and passive (“$+$”) stress by

[TABLE]

[TABLE]

In this model the earth pressure coefficients $K_{x}$, $K_{y}$ are assumed to be functions of the velocity gradient.

The terms $s_{x}$, $s_{y}$ are the net driving accelerations in the $x$, $y$ directions, respectively,

[TABLE]

[TABLE]

where $|\bm{u}|=\sqrt{u^{2}+v^{2}}$, $\kappa=-\frac{\partial\zeta}{\partial x}$ is the local curvature of the reference surface, and $\lambda\kappa$ is the local stretching of the curvature.

The two dimensional Savage-Hutter equations (1)-(3) can be collected in a general vector form as

[TABLE]

where $\bm{U}$ denotes the vector of conservative variables, $\bm{\mathcal{F}}=(\bm{F},\bm{G})$ represent the physical fluxes in the $x$ and $y$ directions, respectively, and $\bm{S}$ is the source term. They are

[TABLE]

The flux for (4) and (5) along a direction $\bm{n}=(n_{x},n_{y})$ is $\bm{\mathcal{F}}\cdot\bm{n}=\bm{F}n_{x}+\bm{G}n_{y}$, while $\bm{n}$ is a unit vector. The Jacobian of the flux along $\bm{n}$ is given by

[TABLE]

where $n_{x}$ and $n_{y}$ are the components of the unit vector in the $x$ and $y$ directions, respectively. The three eigenvalues of the flux Jacobian are

[TABLE]

where $u_{\bm{n}}=\bm{u}\cdot\bm{n}$.

For laboratory avalanches, the ranges of $\delta$ and $\phi$ are usually made $\beta_{x}$ and $\beta_{y}$ greater than zero [11, 23], thus, the eigenvalues of the flux Jacobian are all real and are distinct if the avalanche depth $h>0$. Therefore the Savage-Hutter system given by (4) and (5) is hyperbolic if the gradient of velocity components are not zeros.

Clearly, the Savage-Hutter equations have a very similar mathematical structure to the shallow water equations of hydrodynamics at a first glance. But as a matter of fact, the derivation of Savage-Hutter equations are not the same as those in the shallow water approximation although they both start from the incompressible Navier-Stokes equation. Meanwhile, on account of the jump in the earth pressure coefficients $K_{x_{act/pass}}$ and $K_{y_{act/pass}}$, the complex source term $s_{x}$ and $s_{y}$, the free surface at the front and rear margins, it can be quite complicated to develop an appropriate numerical method to solve the Savage-Hutter equations.

It is a direct observation that the Savage-Hutter equations is not Galinean invariant. Precisely, it is not rotational invariant if $\beta_{x}\neq\beta_{y}$. It is well-known that the shallow water equations is rotational invariant. This is an essential difference between the shallow water equations and the Savage-Hutter equations. Let us point out this fact by

Theorem 1**.**

For any unit vector $\bm{n}=(\cos{\theta},\sin{\theta})$, $\theta\neq 0$, and all vectors $\bm{U}$, the following equality

[TABLE]

holds if and only if $\beta_{x}=\beta_{y}$, where the rotation matrix $\bm{T}$ is as

[TABLE]

Proof.

First we calculate $\bm{TU}$. The result is

[TABLE]

Next we compute $\bm{F}(\bm{TU})$ and obtain

[TABLE]

Then we apply the inverse rotation $\bm{T}^{-1}$ to $\bm{F}(\bm{TU})$ and get

[TABLE]

For the left hand of Eq. (6), we have

[TABLE]

Thus, if $\beta_{x}\neq\beta_{y}$, (7) is not equal to (8). ∎

If the system is rotational invariant, one may solve a 1D Riemann problem with $(h,hu_{\bm{n}})$, where $u_{\bm{n}}=\bm{u}\cdot\bm{n}$, as variable and get a 1D numerical flux $(f_{h},f_{m})$, then the numerical flux along $\bm{n}$ for 2D problem is given as $(f_{h},f_{m}\bm{n})$. This procedure is very popular while it only works in case that $\bm{\mathcal{F}}\cdot\bm{n}=\bm{T}^{-1}\bm{F}(\bm{TU})$ holds. The theorem tells us that this equality holds only when $\theta=0$ or $\beta_{x}=\beta_{y}$. In a finite volume method to be investigated, $\bm{n}$ will be set as the unit normal of the cell interfaces, that it can not always vanish $\theta$ in unstructured grids. Therefore, if $\beta_{x}\neq\beta_{y}$, any approximate Riemann solver based on such a procedure is out of the scope of our candidate numerical flux.

3. Numerical Method

Let us discretize the Savage-Hutter equations (4) and (5) on the Delaunay triangle grids using a finite volume Godunov-type approach. We will adopt the modified HLL approximate Riemann Solver to calculate numerical fluxes across cell interfaces. The MUSCL-Hancock method and ENO-type reconstruction technique are adopted to achieve second-order accuracy in space and time.

3.1. Finite volume method

Before introducing the cell-centered finite volume method, the entire spatial domain $\Omega$ is subdivided into $N$ triangular cells $\tau_{i},i=1,2,\cdots N$. In each cell $\tau_{i}$, the conservative variables of Savage-Hutter equations are discretized as piecewise linear functions as

[TABLE]

where $\bm{U}_{i}^{n}=\left(h_{i}^{n},(hu)_{i}^{n},(hv)_{i}^{n}\right)^{T}$ is the cell average value, $\nabla\bm{U}_{i}^{n}$ is the slope on $\tau_{i}$ and $\bm{x}_{i}$ is the barycenter of $\tau_{i}$, the superscript $n$ indicates at time $t^{n}$.

To introduction the finite volume scheme, Eq. (4) is integrated in a triangle cell $\tau_{i}$

[TABLE]

Applying Green’s formulation, Eq.(9) becomes

[TABLE]

where $\partial\tau_{i}$ denotes the boundary of the cell $\tau_{i}$; $\bm{n}$ is the unit outward vector normal to the boundary; $ds$ is the arc elements. The integrand $\bm{\mathcal{F}}\cdot\bm{n}$ is the outward normal flux vector in which $\bm{\mathcal{F}}=\left[\bm{F},\bm{G}\right]$. Define the cell average:

[TABLE]

where $|\tau_{i}|$ is the area of the cell $\tau_{i}$. The method achieves second-order accuracy by using the MUSCL-Hancock method which consists of predictor and corrector steps. The procedure of the predictor-corrector Godunov-type method goes as follows:

Step 1: Predictor step

[TABLE]

where $\partial\tau_{i,j}$ is the $j$-th edge of $\tau_{i}$ with the unit outer normal as $\bm{n}_{j}$, and $\bm{S}_{i}$ is source term discretized by a centered scheme. The flux vector $\bm{\mathcal{F}}(\bm{U}_{in}^{n})$ is calculated at each cell face $\partial\tau_{i,j}$ after ENO-type piecewise linear reconstruction. $\Delta t$ is the time step.

Step 2: Corrector step

[TABLE]

where the numerical flux vector $\bm{\mathcal{F}}^{*}(\bm{U}_{in}^{n+1/2},\bm{U}_{out}^{n+1/2})\cdot\bm{n}$ is evaluated based on an approximate Riemann solver at each quadrature point on the cell boundary. In the paper, we use HLL approximate Riemann solver owe to its simplicity and stable performance on wet/day interface.

3.2. HLL approximate Riemann solver

For a lot of hyperbolic systems, it uses the approximate Riemann solvers to obtain the approximate solutions in most practical situations. Since Savage-Hutter equations and shallow water equations are very similar, we can directly extend approximate Riemann solvers of shallow water equations to Savage-Hutter equations. The HLL approximate Riemann solver can offer a simple way of dealing with dry bed situaitons and the determination of the wet/dry front velocities [5]. Similar to shallow water equations, the numerical interface flux of HLL for Savage-Hutter equations is computed as follows

[TABLE]

where $\bm{U}_{L}=(h_{L},(hu)_{L},(hv)_{L})$ and $\bm{U}_{R}=(h_{R},(hu)_{R},(hv)_{R})$ are the left and right Riemann states for a local Riemann problem, respectively. $s_{L}$ and $s_{R}$ are estimate of the speeds of the left and right waves, respectively. The key to the success of the HLL approach is the availiability of estimates for the wave speeds $s_{L}$ and $s_{R}$. For the Savage-Hutter equations, there are also dry/wet bed cases. Thus, in order to take the dry bed into account, similar to shallow water equations, the left and right wave speeds are expressed as

[TABLE]

[TABLE]

where $\bm{u}_{L}=(u_{L},v_{L})^{T}$, $\bm{u}_{R}=(u_{R},v_{R})^{T}$, $c_{L}=(\beta_{x}n_{x}^{2}+\beta_{y}n_{y}^{2})_{L}$, $c_{R}=(\beta_{x}n_{x}^{2}+\beta_{y}n_{y}^{2})_{R}$ and

[TABLE]

3.3. Linear reconstruction

In order to get second order accuracy in spatial, a simple ENO-type piecewise linear reconstruction was applied to suppress the numerical oscillations near steep gradients or discontinuites. This reconstrction method had been used by Deng et. al. for shallow water equations [3].

This part is mainly taken from [3]. For a variable $\psi=h,hu$ or $hv$, assuming its cell average value $\psi_{i}$ on $\tau_{i}$ is given (see Fig. 2), then we can get the gradients of $\psi$ for the three triangles $\triangle AiB$, $\triangle BiC$ and $\triangle CiA$ are referred to as $(\nabla\psi)_{AiB}$, $(\nabla\psi)_{BiC}$ and $(\nabla\psi)_{CiA}$, and choose the one with the minimal $l^{2}$ norm as $\nabla\psi_{i}$, that is,

[TABLE]

Acoording to the idea of minmod slope limiter, we let $\nabla\psi_{i}=0$ when $\psi_{i}\geq\max{\{\psi_{A},\psi_{B},\psi_{C}\}}$ or $\psi_{i}\leq\min{\{\psi_{A},\psi_{B},\psi_{C}\}}$. Meanwhile, in the process of reconstructing $h$, we must consider the physical criteria, that is, the avalanche depth $h$ must be non-negative at each quadrature points, otherwise $\nabla h_{i}$ is set to be zero.

3.4. Time stepping

The time-marching formula is explicit, and the time step length $\Delta t$ is determined by the Courant-Friedrichs-Lewy(CFL) condition for its stability. For triangular grid and implemention of CFL condition, the time step length is usually expressed as

[TABLE]

where $\mathrm{C_{r}}$ is the Courant number specified in the range $0<\mathrm{C_{r}}\leq 1$, and $\mathrm{C_{r}}=0.5$ is adopted in our simulations; $|\bm{\beta}|=\sqrt{\beta_{x}^{2}+\beta_{y}^{2}}$, $N$ is the total number of cells, and $\Delta x_{i}$ is the size of $\tau_{i}$. For triangular grid, $\Delta x_{i}$ is usually calculated as the minimum barycenter-to-barycenter distance between $\tau_{i}$ and its adjacent cells.

4. Numerical Examples

In this section, three numerical examples were performed to verify the proposed scheme which implemented using C++ programming language based on the adaptive finite element package AFEPack[16]. The initial triangular grids were generated by easymesh. Because the posteriori error estimator for $h$-adaptive method on unstructured grid was not the main goal of this work, we just applied the following local error indicator which advised by [3]

[TABLE]

where

[TABLE]

in which $\tau$ is a cell in the mesh and $|\tau|$ is its area, $\left[[\cdot\right]]$ denotes the jump of a variable across $\partial\tau$, and $I_{h}\psi$ is the piecewise linear numerical solution for $\psi=h,hu,hv$, respectively.

4.1. Dam break test cases with exact solution

For this dam break problem, it was originally a one dimensional problem with analytical solutions involving a constant bed slope and a Coulomb frictional stress, thus, its governing equations can reduced as [14][30]

[TABLE]

where $g$ is the gravitational constant, the analytical solution of a granular dam break problem is given [30]:

[TABLE]

where $c_{0}=\sqrt{gh_{0}\cos{\zeta}}$ and

[TABLE]

[TABLE]

In order to check the performance of the numerical schemes proposed in this work, we modified it to an equivalent two-dimensional problem as Ref.[30]. The computational domain is $[-12.8,12.8]\times[-1.6,1.6]$, the initial conditions are

[TABLE]

and $\zeta=40^{\circ},\delta=24.5^{\circ}$.

Fig. 3 plots the compute avalanche depth in comparison with the exact solutions at $t=0.1,0.2,0.3,0.4$ and $0.5$. It can be seen that our numerical solutions are in agreement with exact solutions very well, which also verifies the robustness and validity of our numerical algorithm.

4.2. Avalanche flows slide down an inclined plane and merge

continuously into a horizontal plane

In this section, we simulate an avalanche of finite granular mass sliding down an inclined plane and merging contiuously into a horizontal plane, which was considered in [29] to verify their model’s ability.

In order to compare with the results from Ref. [29][30], all computational parameters are set to the same as they are. That is, the computational domain is $[0,30]\times[-7,7]$, and $\epsilon=\lambda=1$ in our simulations. The inclined section lies $x\in[0,17.5]$ and the horizontal region lies $x\in[21.5,30]$ with a smooth transition zone lies $x\in[17.5,21.5]$. The inclination angle is given by

[TABLE]

where $\zeta_{0}=35^{\circ}$ and $\delta=\phi=30^{\circ}$. The granular mass is suddenly released at $t=0$ from the hemispherical shell with an initial radius of $r_{0}=1.85$ in dimensionless length units. The center of the cap is initially located at $(x_{0},y_{0})=(4,0)$.

Fig. 4 shows the avalanche depth contours of the fluid at eight time slices as the avalanche slides on the inclined plane into the horizontal run-out zone. From the Fig. 4 it can be seen that, the avalanche starts to flow along the $x$ and $y$ direction due to gravity, and it flows faster along the downslope direction obviously. At $t=9$ the part of granular reaches the horizontal run-out zone and begins to deposite because of the basal friction is sufficient large, and the tail is still doing the acceleration movement. Here, some small numerical oscillations are visible in the solutions which also found in Ref.[29][30] with orthogonal quadrilateral mesh. At $t=12$ a shock wave develops round the end of the transition zone at $x=21.5$ and a surge wave is created. During the $t=12$ to $t=24$, most of graunlar mass are accumulated at the end of the transition zone and the front-end of the horizontal zone, meanwhile, a obvious backward surge can be found. At $t=24$ the final deposition of the avalanche is nearly formed and an expansion fan is observed.

On the whole, our numerical resluts are identical with these available from Ref. [29][30]. Fig. 5 presents the mesh adaptive moving for the avalanche granular flows down the inclined plane and merges contiuously into a horizontal plane at $t=6,12,18,24$. It can be seen that the triangle mesh is refined at the free surface where gradients varied significantly.

4.3. Granular flow pasts obstacle

The research of granular avalanches around obstacle is of partiacular interest because some important pheomena such as shock waves, expansion fans and granular vacuum can be observed [2]. Meanwhile, setting obstacles is an effective means to control the dynamic process of granular avalanche flow. In this section, we present a simulation example of finite granular mass flows down an inclined plane and merges inito a horizontal run-out zone with a conical obstacle located on the inclined plane. In addition, in order to compared the effects of with and without obstacle, all compuational parameters (such as computational domain, the inclination angle $\zeta(x)$, $\phi$, $\delta$ ) are the same as the above example. A conical obstacle is incorporated into this model and the centre of obstacle is located at $(13,0)$ with dimensionless radius of bottom-circle of 1 and dimsionless heights $H=1$.

Fig. 6 shows a series of snap-shots of the thickness contours as the avalanche slides down the inclined plane which has a conical obstacle in the avalanche track. It can be seen that the flow of granular is the same as that without obstacle before it reaches the obstacle (eg. for $t\leq 3$). At $t=6$ the granular avalanche has reache the obstacle, some of them climb the obstacle, and the other divide into two parts and flow around the obstacle downwards. At the moment, a stationary zone appears in the front of the obstacle, where the avalance velocity becomes zero, but it is followed by a great increase in the avalanche thickness. Meanwhile, behind the obstacle, a so-called granular vacuum is formed which can protect the zone directly behind it from disasters. At $t=9$, the granular mass continues to accelerate until it reaches the horizontal run-out zone where the basal friction is sufficient large, which results in the front becomes to rest, but the part of the tail accelerates further. At $t=12$ the shcoks waves fromed on the horizontal run-out zone during the front of granular mass deposits and the tail of it still slides accelerated. From $t=15$ to $t=18$, most of the granualr mass has slided down either side of the obstacle and on the transition zone two expansion fans gradually formed. For $t=21,24$ two separate and symmetrical deposition mass of granular are formed at the end of transition zone and the front of horizontal zone. From $t=6$ to $t=24$ there always has the flow-free region directly behind the obstacle, which can be regarded as the protected area in the practial application of avalanche protection. Thus, the obstacle does change the path of granular avalanche flow movement and affects its final deposition pattern. The numerical results show that our numerical model has the ability of capturing key qualitative features, such as shocks wave and flow-free region. As in the previous example, Fig. 7 also shows the adaptive mesh moving for the granular avalanche flows down with a obstacle at $t=6,12,18,24$.

5. Conclusions

This paper is an attempt to solve the Savage-Hutter equations on unstructred grids. Keeping in mind that the model is not rotational invariant, it is more than encouraging for our numerical scheme to attain numerical results with good agreement to the reference solutions. It is indicated that the method is applicable to problems as granular avalanche flows. The underlying reason why the scheme works smoothly is still under investigation.

Acknowledgement

This work is financially supported by the National Natural Science Foundation of China(No. 11826207 and No. 11826208).