A fifth order WENO scheme for numerical simulation of shallow granular two-phase flow model

Abstract

In this article, a weighted essentially non-oscillatory (WENO) scheme is implemented to simulate two-phase shallow granular flow (TPSF) model. The flow is assumed to be incompressible and it is regarded as shallow layer of granular and liquid material. The mathematical model consists of two phases, that is, solid and liquid. Each phase has its continuity and momentum equation. The presence of the equations are coupled together involving the derivatives of unknowns which make it more challenging to solve. An efficient numerical technique is needed to tackle the numerical complexities. Our main intrigue is the numerical approximation of the above-mentioned solid-liquid model. The weighted essentially non-oscillatory (WENO) scheme of order 5 is utilized to handle the shock waves and contact discontinuities appear in the solution. The results are compared with the results already available in the literature by conservation element and solution element (CESE) scheme. It is observed the WENO scheme produces less errors as compared to CESE scheme and also effectively handle the shocks.

Keywords

Shallow granular two-phase flow model WENO scheme CESE scheme non-conservative systems

Introduction

In this article, we investigate depth-averaged two-phase shallow granular flow (TPSF) model for gravity-driven mixtures of solid and liquid. The motivation behind studying two-phase flows is because of its wide range applications to debris flows and avalanches.¹ In most of the cases, the mathematical equations of these flows form hyperbolic conservation laws. The most challenging task for the scientists is the approximation of such flows. Noteworthy endeavors have been made in recent years for modeling and simulation of two-phase flows.

Savage and Hutter^2,3 commenced the work on one-dimensional single-phase flow model which was proposed to study avalanches. Later, they extended it to two-dimensions⁴ and several authors generalized the results on more complex basal topography.^5–10 Through depth-averaged models, incredible progress has been made in enhancing numerical and scientific evidence of geophysical flows. As of late, numerous higher order numerical schemes implementing these depth-averaged models have been able to solve with significant success.^11–15 The two-phase shallow granular flow model considered in this article is a variant of Pitman and Le¹⁶ fluid model, comprising of shallow layer of liquid and solid phases. In literature, the impacts of the liquid phase^17–19 are not taken into consideration in a number of models used to approximate true avalanches. Though it plays a crucial role in the mobility of flow.²⁰

The equations of two-phase shallow granular flow model form a system of hyperbolic conservation laws. This model have been studied by numerous scientists while assuming sufficiently small differences between phase velocities.^{14,15,21–23} When analyzing geophysical flows, the assumption of hyperbolicity is reasonable because solid and liquid are quite often quickly brought into kinematic equilibrium by inter-phase interactions.²⁰ Pelenti et al.^14,15 approximated the two-phase granular model by using a finite volume Godunov-type technique based on the Riemann solver (Roe-type). Unfortunately, the Roe-type scheme produces negative values of physical variables like height and solid volume fraction. Later, Riemann solver based finite volume scheme (FVS) was implemented to solve the two-phase shallow granular flow model.²¹ In Yan et al.,²⁴ the authors numerically investigated the Eulerian-Eulerian and Eulerian-Lagrangian two-phase flow models. The kinetic flux-vector splitting (KFVS) and CE/SE schemes were extended to solve two-phase shallow granular flow.^22,23 In Zia and Qamar,²³ the schemes preserves the positivity of height only for $γ = \frac{1}{2}$ . Therefore, our motivation was to extend a scheme which can handle different flow conditions for all values of $γ < 1$ (ratio of densities of two phases).

A well-balanced numerical scheme can handle small perturbations in the steady state solutions effectively in comparison with ordinary numerical scheme (not well-balanced). Therefore well-balanced schemes are more appropriate to compute turbulent fluctuations mainly in steady flows, whereas schemes which does not possess the property of being well-balanced can only handle perturbations at the level of truncation error with the given grid. For the depth averaged models, a lot of work in literature have been done on well-balanced schemes.^{14,15,21,25–28}

In this study, a fifth-order finite volume WENO numerical technique²⁹ is designed to solve the shallow two-phase flow model. The proposed numerical scheme is more simple and fundamentally distinctive than the computational methods that were devised to solve the mentioned model in the articles.^14,15 Initially in 1987, Harten and Osher³⁰ introduced the finite volume essentially non-oscillatory (ENO) numerical scheme. This scheme obtained arbitrarily high order accuracy in smooth regions, resolved the steep gradients efficiently and did not produce spurious oscillations in the vicinity of sharp gradients. Later in 1988, Shu and Osher³¹ developed the finite difference ENO scheme with total variation diminishing Runge-Kutta (TVD-RK) time discritization procedure. Subsequently, Liu et al.³² proposed the improved version of finite volume ENO numerical scheme, namely third order finite volume weighted essentially non-oscillatory (WENO) numerical scheme. Soon after in 1996, Jiang and Shu³³ introduced the fifth order WENO numerical scheme. After that, the WENO numerical approach has been revised, extended, and improved for several fields of science and engineering, for detail see References.^34–44 The key idea in developing the WENO numerical scheme is used a convex combination of interpolations or reconstructions from different candidate stencils, rather than using only one of them as in the case of ENO. The combination coefficients, also called nonlinear weights, are chosen in such a way that if all candidate stencils contain smooth solution then these nonlinear weights should be selected very close to the linear weights, which could provide the highest possible order of accuracy from the combined stencil. Further, if a candidate stencil contains a steep gradient, while at least one other candidate stencil that does not contain steep gradient, then the candidate stencil that contains a steep gradient should be assigned a very small weight, thus its influence toward the approximation is negligible, for more detail see Shu et al.^33–35 Hence, the suggested numerical approach effectively resolves sharp discontinuities while maintaining high order accuracy in smooth areas.

This article is structured as follows. Section 2 presents 1D shallow granular two-phase flow model. Section 3 is devoted for the derivation of WENO scheme. The numerical test cases are carried out in section 4. Finally, conclusions are drawn in section 5.

Shallow granular two-phase flow model

A shallow layer of liquid-solid mixture over an even surface is considered. Let $ρ_{f}$ and $ρ_{s}$ be the densities of liquid and solid phases of an incompressible flow respectively. The height of the flow be $h$ and $ϕ$ be the solid volume fraction. The other variables are defined as

h_{s}^{x} = ϕ h and h_{f}^{x} = (1 - ϕ) h .

(1)

The velocities in $x$ -direction of both solid and fluid phases are denoted by $u_{s}^{x}, u_{f}^{x}$ , respectively. Phase momenta are given by $m_{s}^{x} = h_{s}^{x} u_{s}^{x}$ and $m_{f}^{x} = h_{f}^{x} u_{f}^{x}$ . The continuity and momentum equations for the two-phase shallow system are^21,22

\partial_{t} h_{s}^{x} + \partial_{x} m_{s}^{x} = 0,

(2)

\begin{matrix} \partial_{t} m_{s}^{x} + \partial_{x} (\frac{{(m_{s}^{x})}^{2}}{h_{s}^{x}} + \frac{g}{2} {(h_{s}^{x})}^{2} + \frac{g}{2} (1 - γ) h_{s}^{x} h_{f}^{x}) \\ + γ {gh}_{s}^{x} \partial_{x} h_{f}^{x} = - {gh}_{s}^{x} \partial_{x} b + γ F^{D}, \end{matrix}

(3)

\partial_{t} h_{f}^{x} + \partial_{x} m_{f}^{x} = 0,

(4)

\partial_{t} m_{f}^{x} + \partial_{x} (\frac{{(m_{f}^{x})}^{2}}{h_{s}^{x}} + \frac{g}{2} {(h_{f}^{x})}^{2}) + {gh}_{f}^{x} \partial_{x} h_{s} = - {gh}_{f}^{x} \partial_{x} b - F^{D} .

(5)

Here, $b \overset{def}{=} b (x)$ is basal topography, $g$ is the gravity constant, $F^{D}$ is representing inter-phase drag forces which are given by, $F^{D} = D (h_{s}^{x} + h_{f}^{x}) (u_{f}^{x} - u_{s}^{x})$ , where $D$ is a drag function depending generally upon $ϕ$ , $| u_{f}^{x} - u_{s}^{x} |$ and a set of physical parameters $σ$ (e.g. particle diameter, specific densities) and $γ = ρ_{f} / ρ_{s} < 1$ . If we take $D = 0$ , then the equations above are somewhat similar to the two-layer shallow flow models.^45–47 In this article, we do not consider the drag forces. In compact form the system of equations are expressed as²¹

\frac{\partial W}{\partial t} + \frac{\partial F (W)}{\partial x} = S (W),

(6)

where,

\begin{matrix} W = (\begin{matrix} h_{s}^{x} \\ m_{s}^{x} \\ h_{f}^{x} \\ m_{f}^{x} \end{matrix}), F (W) = (\begin{matrix} m_{s}^{x} \\ \frac{{(m_{s}^{x})}^{2}}{h_{s}^{x}} + \frac{g}{2} {(h_{s}^{x})}^{2} + \frac{g}{2} (1 - γ) h_{s}^{x} h_{f}^{x} \\ m_{f}^{x} \\ \frac{{(m_{f}^{x})}^{2}}{h_{s}^{x}} + \frac{g}{2} {(h_{f}^{x})}^{2} \end{matrix}), S (W) = (\begin{matrix} 0 \\ - {gh}_{s}^{x} \partial_{x} b \\ 0 \\ - {gh}_{f}^{x} \partial_{x} b \end{matrix}), \end{matrix}

(7)

The continuity equations for $h_{s}^{x}$ and $h_{f}^{x}$ are conservative, whereas $m_{s}^{x}$ and $m_{f}^{x}$ representing momentum equations for both solid and liquid phases respectively. The system can also be written in following quasi-linear form

\partial_{t} W + A (W) \partial_{x} (W) = η^{b} (W),

(8)

where, $η^{b} (W) = {(0, - {gh}_{s}^{x} \partial_{x} b, 0, - {gh}_{f}^{x} \partial_{x} b)}^{T}$ and $A (W)$ is given by

A (W) = (\begin{matrix} 0 & 1 & 0 & 0 \\ - {(\frac{q_{2}}{q_{1}})}^{2} + g q_{1} + g \frac{1 - γ}{2} q_{3} & 2 \frac{q_{2}}{q_{1}} & g \frac{1 + γ}{2} q_{1} & 0 \\ 0 & 0 & 0 & 1 \\ g q_{3} & 0 & - {(\frac{q_{4}}{q_{3}})}^{2} + g q_{3} & 2 \frac{q_{4}}{q_{3}} \end{matrix}) .

(9)

The conditions on the variables for this model are:

D = {W \in R^{4}; h_{s}^{x}, h_{f}^{x} \geq 0 and u_{s}^{x}, u_{f}^{x} \in R},

(10)

or, with $h = h_{s}^{x} + h_{f}^{x}$ and $ϕ = \frac{h_{s}^{x}}{h_{s}^{x} + h_{f}^{x}}$ , the conditions listed above can be written as follows:,

D = {W \in R^{4}; h \geq 0, 0 \leq ϕ \leq 1, and u_{s}^{x}, u_{f}^{x} \in R} .

(11)

Introduction to fifth order WENO Scheme

The development of a finite volume WENO numerical technique for numerical approximation of two-phase flow model is presented in this section. Because the model under consideration can be expressed in a concise manner as

\partial_{t} W + \partial_{x} F (W) = S (W, x), x \in Ω, t > 0

(12)

as well as dividing the space $Ω$ into cells $C_{i} = [x_{i - \frac{1}{2}}, x_{i + \frac{1}{2}}], i = 1, 2, 3 . . ., N .$ The center and length of the $i - th$ cell is indicated by in this case by $x_{i} = \frac{(x_{i - \frac{1}{2}} + x_{i + \frac{1}{2}})}{2}$ and $Δ x_{i}$ respectively. Upon integration of equation (12) over $C_{i}$ , we get

\begin{matrix} \frac{d}{dt} W_{i} (t) + \frac{1}{Δ x_{i}} (F (W (x_{i + \frac{1}{2}}, t) - F (W (x_{i - \frac{1}{2}}, t)) \end{matrix}) = \frac{1}{Δ x_{i}} \int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} S (W, x) dx,

(13)

where $W_{i} (t) = \frac{1}{Δ x_{i}} \int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} W (x, t) dx$ . Approximation of equation (13) leads to,

\frac{d}{d t} W_{i} (t) + \frac{1}{Δ x_{i}} ({\hat{F}}_{i + \frac{1}{2}} - {\hat{F}}_{i - \frac{1}{2}}) = \frac{1}{Δ x_{i}} \int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} S (W, x) d x .

(14)

Where, monotone numerical flux is denoted by ${\hat{F}}_{i + \frac{1}{2}} = Ϝ (W_{i + \frac{1}{2}}^{-}, W_{i + \frac{1}{2}}^{+})$ and $W_{i + \frac{1}{2}}^{\pm}$ are pointwise approximations to $W (x_{i + \frac{1}{2}}, t)$ . As a monotone numerical flux, we shall utilize the Lax-Friedrichs flux (LFF), which is defined as follows

\begin{array}{l} Ϝ (W_{i + \frac{1}{2}}^{-}, W_{i + \frac{1}{2}}^{+}) = \frac{1}{2} \\ (F (W_{i + \frac{1}{2}}^{-}) + F (W_{i + \frac{1}{2}}^{+}) - ϑ (W_{i + \frac{1}{2}}^{+} - W_{i + \frac{1}{2}}^{-})), \end{array}

(15)

where $ϑ = max_{W} | F' (W) |$ . Cell averages $W (x_{i}, t)$ are approximated by the computational variables $W_{i} (t)$ . By using WENO reconstruction, the $W_{i + \frac{1}{2}}^{-}$ are computed from the surrounding cell-average values B. By using WENO reconstruction, the $W_{i + \frac{1}{2}}^{-}$ and $W_{i + \frac{1}{2}}^{+}$ are derived using the nearby cell-average values $W_{i}$ . The following relationships yield the point-wise reconstructed values $W_{i + \frac{1}{2}}^{+}$ and $W_{i + \frac{1}{2}}^{-}$ :

W_{i + \frac{1}{2}}^{+} = ω_{1} {\hat{W}}_{i + \frac{1}{2}}^{1} + ω_{2} {\hat{W}}_{i + \frac{1}{2}}^{2} + ω_{3} {\hat{W}}_{i + \frac{1}{2}}^{3},

(16)

W_{i + \frac{1}{2}}^{-} = {\tilde{ω}}_{1} {\tilde{W}}_{i + \frac{1}{2}}^{1} + {\tilde{ω}}_{2} {\tilde{W}}_{i + \frac{1}{2}}^{2} + {\tilde{ω}}_{2} {\tilde{W}}_{i + \frac{1}{2}}^{2},

(17)

where ${\hat{W}}_{i + \frac{1}{2}}^{l}$ and ${\tilde{W}}_{i - \frac{1}{2}}^{l}$ , for $l = 1, 2, 3$ are values that have been recreated and described as

{\hat{W}}_{i + \frac{1}{2}}^{1} = \frac{1}{6} (2 W_{i}^{+} + 5 W_{i + 1}^{+} - W_{i + 2}^{+}),

(18)

{\hat{W}}_{i + \frac{1}{2}}^{2} = \frac{1}{6} (- W_{i - 1}^{+} + 5 W_{i}^{+} + 2 W_{i + 1}^{+}),

(19)

{\hat{W}}_{i + \frac{1}{2}}^{3} = \frac{1}{6} (2 W_{i - 2}^{+} - 7 W_{i - 1}^{+} + 11 W_{i}^{+}),

(20)

and

{\tilde{W}}_{i + \frac{1}{2}}^{1} = \frac{1}{6} (11 W_{i + 1}^{-} - 7 W_{i + 2}^{-} + W_{i + 3}^{-}),

(21)

{\tilde{W}}_{i + \frac{1}{2}}^{2} = \frac{1}{6} (W_{i}^{-} + 5 W_{i + 1}^{-} - W_{i + 2}^{-}),

(22)

{\tilde{W}}_{i + \frac{1}{2}}^{3} = \frac{1}{6} (- W_{i - 1}^{-} + 5 W_{i - 1}^{-} + W_{i}^{-}) .

(23)

The nonlinear weights $ω_{l}$ and ${\tilde{ω}}_{l}$ in equations (16) and (17) are defined as

ω_{l} = \frac{ξ_{l}}{\sum_{m = 1}^{3} ξ_{m}} with ξ_{l} = \frac{σ_{l}}{ϵ + B_{l}}, l = 1, 2, 3,

(24)

and

{\tilde{ω}}_{l} = \frac{{\tilde{ξ}}_{l}}{\sum_{m = 1}^{3} {\tilde{ξ}}_{m}} with {\tilde{ξ}}_{l} = \frac{{\tilde{σ}}_{l}}{ε + {\tilde{B}}_{l}}, l = 1, 2, 3,

(25)

To avoid the denominator becoming zero, $ε$ can be any tiny positive number, and we set $ε = 10^{- 7}$ in all approximations.

The linear weights are denoted by $σ_{l}$ and ${\tilde{σ}}_{l}$ , whereas the smoothness indicators are denoted by $B_{l}$ and ${\tilde{B}}_{l}$ . $σ_{l}$ and ${\tilde{σ}}_{l}$ ’s values are as follows:

σ_{l} = {\begin{matrix} \frac{3}{10} if l = 1, \\ \frac{6}{10} if l = 2, \\ \frac{1}{10} if l = 3, \end{matrix} and {\tilde{σ}}_{l} = {\begin{matrix} \frac{1}{10} if l = 1, \\ \frac{6}{10} if l = 2, \\ \frac{3}{10} if l = 3 . \end{matrix}

(26)

$B_{l}$ and ${\tilde{B}}_{l}$ are smoothness indicators and have the following values:

B_{l} = {\begin{matrix} \frac{13}{12} {(W_{i}^{+} - 2 W_{i + 1}^{+} + W_{i + 2}^{+})}^{2} \\ + \frac{1}{4} {(W_{i}^{+} - 4 W_{i + 1}^{+} + W_{i + 2}^{+})}^{2} if l = 1, \\ \frac{13}{12} {(W_{i - 1}^{+} - 2 W_{i}^{+} + W_{i + 1}^{+})}^{2} \\ + \frac{1}{4} {(W_{i - 1}^{+} - W_{i + 1}^{+})}^{2} if l = 2, \\ \frac{13}{12} {(W_{i - 2}^{+} - 2 W_{i - 1}^{+} + W_{i}^{+})}^{2} \\ + \frac{1}{4} {(W_{i - 2}^{+} - 4 W_{i - 1}^{+} + W_{i}^{+})}^{2} if l = 3 . \end{matrix}

(27)

and

{\tilde{B}}_{l} = {\begin{matrix} \frac{13}{12} {(W_{i + 1}^{-} - 2 W_{i + 2}^{-} + W_{i + 3}^{-})}^{2} \\ + \frac{1}{4} {(3 W_{i + 1}^{-} - 4 W_{i + 2}^{-} + W_{i + 3}^{-})}^{2} if l = 1, \\ \frac{13}{12} {(W_{i}^{-} - 2 W_{i + 1}^{-} + W_{i + 2}^{-})}^{2} \\ + \frac{1}{4} {(W_{i}^{-} - W_{i + 2}^{-})}^{2} if l = 2, \\ \frac{13}{12} {(W_{i - 1}^{-} - 2 W_{i}^{-} + W_{i + 1}^{-})}^{2} \\ + \frac{1}{4} {(W_{i - 1}^{-} - 4 W_{i}^{-} + 3 W_{i + 1}^{-})}^{2} if l = 3 . \end{matrix}

(28)

The spatial reconstruction process is now finished. On the right-hand side of equation (13), the nonzero terms is discretized after that,

\int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} - {gh}_{s}^{x} \partial_{x} bdx = \int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} (- {gh}_{i}^{x} \partial_{x} (b)) dx

(29)

= - (\bar{g} h_{i}^{x}) ({(b)}_{i + \frac{1}{2}} - {(b)}_{i - \frac{1}{2}})

(30)

The averaged value at cell $C_{i}$ is $\bar{g} h_{i}^{x}$ , while the source fluxes at the cell $C_{i}$ ’s interfaces are $(b)_{i \pm \frac{1}{2}}$ . In addition, for $C_{i}$ ’s right interface, we write as

(b)_{i + \frac{1}{2}} = (b)_{i}^{+} + (b)_{i + 1}^{-}

(31)

The values of cell $C_{i}$ ’s left interface and $\int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} - {gh}_{f}^{x} \partial_{x} b$ ’s dicretization are determined in the same way.

Further, $S_{i}$ denotes the approximation of the expression $\int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} S (W, x) dx$ . Finally, we arrive at the following semi-discrete equation:

\frac{d}{dt} {\bar{W}}_{i} (t) = - \frac{1}{Δ x_{i}} ({\hat{F}}_{i + \frac{1}{2}} - {\hat{F}}_{i - \frac{1}{2}}) + \frac{1}{Δ x_{i}} S_{i},

(32)

alternatively, the equation above may be written as

\frac{d}{dt} {\bar{W}}_{i} (t) = O (W) .

(33)

Now, we use the third-order TVD RK technique²⁹ to solve the system of ordinary differential equation (33) as

\begin{matrix} W^{(1)} = W^{n} + dt O (W^{n}), \\ W^{(2)} = \frac{3}{4} W^{n} + \frac{1}{4} (W^{(1)} + dt O (W^{(1)})), \\ W^{(n + 1)} = \frac{1}{3} W^{n} + \frac{2}{3} (W^{(2)} + dt O (W^{(2)})), \end{matrix}

(34)

$O (W)$ is the spatial operator in this case. The Courant-Friedrichs-Lewy condition denoted by $CFL$ is used to estimate the time step $dt$ using the relation

\frac{{CFL}^{*} dx}{max (| λ_{1} |, | λ_{2} |, | λ_{3} |, | λ_{4} |, | λ_{5} |)} .

Numerical test problems

To validate the proposed WENO scheme several numerical test problems are considered here. In all the experiments, we set $γ < 1$ . For comparison, we use the results of CESE scheme. The codes for simulating these problems are written in C++and graphs are plotted through MATLAB.²²

Problem 1: $L^{1}$ -error of the schemes.

In this test case, we computed the $L^{1}$ -errors in height $h$ and velocity $u$ for WENO and CESE schemes. The findings of $L^{1}$ -errors are provided in the Table 1. A comparison from the table demonstrates that WENO scheme has less error as compared to CESE scheme.

Problem 2: Simple Riemann test problem.

Table 1.

Problem 3: $L^{1}$ -error comparison.

N	$h$		$u$
	WENO	CESE	WENO	CESE
50	0.0145	0.0132	0.2417	0.0341
100	0.0071	0.0064	0.0094	0.0113
200	0.00021	$0.0032$	0.0006	0.0069
400	0.00031	$0.0017$	0.0002	0.0035
800	$5.0321 \times 10^{- 4}$	$7.8010 \times 10^{- 4}$	$0.0008$	0.0015
1600	$2.0627 \times 10^{- 4}$	$3.3049 \times 10^{- 4}$	$4.6361 \times 10^{- 4}$	$8.6179 \times 10^{- 4}$

This problem is also discussed in References.^15,22,23 The interface is located at $x = 0$ initially. The values of $h, ϕ, u_{s}^{x}, u_{f}^{x}$ at right and left of interface are given below

(h, ϕ, u_{s}^{x}, u_{f}^{x})_{L} = (1, 0.8, 0, 0),

(35)

(h, ϕ, u_{s}^{x}, u_{f}^{x})_{R} = (2, 0.4, - 0.9, - 0.1) .

(36)

Here, we take $b (x) = constant$ , $g = 9.81$ and $L = 10$ . The domain $- 5 \leq x \leq 5$ is discritized into $300$ grid points. At $t = 0.5$ the results for $h$ , $h_{s}^{x}$ , $h_{f}^{x}$ , the solid volume fraction $ϕ$ and the phase velocities $u_{s}^{x}$ , $u_{f}^{x}$ are shown in Figure 1. The approximated solution consists of 1 – rarefaction wave, 2 – shock wave, and 3 – rarefaction wave respectively. On comparison with CESE scheme, WENO provides the better resolution of discontinues profiles. Furthermore, our results are in close agreement with those available in literature.^15,22,23

Problem 3: Dam break problem (Rarefaction into vacuum of the fluid constituent).

Figure 1.

Problem 2: simple Riemann problem: on 300 cells, t = 0.5 yielded the following results.

Here, the flow for $h > 0$ is considered on the entire spatial domain and for all times, but characterized by a nearly vacuum region for the liquid phase $(h_{f} = 0)$ . The Riemann data for this problem is:

(h, ϕ, u_{s}^{x}, u_{f}^{x})_{L} = (1, 0.8, 0, 0), if x \leq 0,

(37)

(h, ϕ, u_{s}^{x}, u_{f}^{x})_{R} = (1, 1, 0, 0), if x > 0 .

(38)

There is no liquid on the right side of the interface $x = 0$ , that is, vacuum zone. The outflow boundary conditions are taken here with $g = 9.81$ . The domain is $[- 5, 5]$ is divided into $300$ cells. To prevent division by zero, we bring on the right hand side $ϕ = 1 - 10^{- 6}$ . Figure 2 shows the outcomes for this problem. We can observe that the solution consists of two-rarefaction waves across which $h_{f} = h (1 - ϕ)$ vanishes. In addition, four more rarefaction waves and one shock wave are observed in a pure solid area where $(h_{f}^{x} = 0)$ is observed. One can also observe that results are free of oscillations.

Problem 4: Spreading of a granular mass

Figure 2.

Problem 3: dam break problem. Results at $t = 1$ .

In this problem, we consider the case of spreading of granular mass on a horizontal surface. The velocities of both the phases are taken as zero initially, that is $(u_{s}^{x} = u_{f}^{x} = 0)$ . The other Riemann data for this problem are,

h (x) = {\begin{matrix} 1, & - 1 \leq x \leq 1; \\ 0, & otherwise . \end{matrix} and ϕ (x) = (3 + 4 e^{- x^{2}}) / 10 .

(39)

Here, $b (x) = constant$ and $g = 1.0$ . In this test case we take $h = 0$ below the tolerance $\in = 10^{- 5}$ . The domain $- 10 \leq x \leq 10$ is divided into $300$ cells. Though inter-phase drag forces are not considered here, still $| u_{s}^{x} - u_{f}^{x} | \to 0$ as $h \to 0$ . The results are depicted in the Figure 3. The hyperbolic conditions are maintained as $h \to 0$ . Also a noticeable change in of $ϕ$ (solid volume fraction) can also be observed.

Problem 5: Propagation at interface.

Figure 3.

Problem 4: spreading of granular mass.

The purpose of considering this test problem is to check the effectiveness of our suggested WENO scheme in capturing the discontinuous profiles over time. This is a very hard test problem for any numerical scheme. Several authors used this test case to validate numerical schemes.^22,28,48 The initial data are given as

(h_{s}^{x}, h_{f}^{x}, u_{s}^{x}, u_{f}^{x})_{L} = (0.5, 0.5, 2.5, 2.5), if x \leq 0.5,

(40)

(h_{s}^{x}, h_{f}^{x}, u_{s}^{x}, u_{f}^{x})_{R} = (0.55, 0.45, 2.5, 2.5), if x > 0.5 .

(41)

Here, $g = 9.81$ and $γ = 0.98$ . The initial interface is located at $x = 0.5$ . The domain $0 \leq x \leq 1$ is discretized into $500$ grid points. The solution involves the interface evolution between the two phases, each phase velocity and free surface displacement see Figure 4. The WENO scheme effectively handled the sharp discontinuous profiles.

Problem 6: Dry bed generation.

Figure 4.

Problem 5: interface evolution outcomes (top left), free surface displacement (top right), upper layer velocity (bottom left) and lower layer velocity (bottom right).

This test problem is related to the formation of a dry bed zone. The left and right state data for this problem is,

(h, ϕ, u_{s}^{x}, u_{f}^{x})_{L} = (0.1, 0.4, - 3.0, - 3.0), if x \leq 0,

(42)

(h, ϕ, u_{s}^{x}, u_{f}^{x})_{R} = (0.1, 0.7, 3.0, 3.0), if x > 0 .

(43)

The initial discontinuity is at $x = 0$ and $b (x) = constant$ . The gravitational constant $g = 9.81$ and $γ = 0.2$ . The computational domain $[- 5, 5]$ is divided into 100 points.

The results are shown in Figure 5. The solution consist of two oppositely moving rarefaction waves generating a dry bed region in between. The results by WENO scheme are compared with the CESE scheme. The WENO scheme gives better resolution of the shocks.

Figure 5.

Problem 6: dry bed generation.

Conclusions

An incompressible two-phase shallow granular flow model is numerically investigated by fifth order WENO scheme. The proposed scheme is capable of capturing sharp discontinuities and accurately handle the dry bed zones for all values of $γ < 1$ . The test problems considered here show the ability of the suggested technique to tackle a variety of flow conditions involving vacuum formation and dry beds. The results by WENO scheme are compared with those obtained from CESE scheme. A good agreement was observed between the two schemes. However, WENO scheme produces less error than CESE scheme and accurately resolve the sharp discontinuous profiles.

Footnotes

Handling Editor: Chenhui Liang

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iDs

Omar Rabbani

Saqib Zia

Asad Rehman

References

Toro

EF.

Shock-capturing methods for free-surface shallow flows. New York, NY: John Wiley, 2001.

Savage

Hutter

The motion of a finite mass of granular material down a rough incline. J Fluid Mech 1989; 199: 177–215.

Savage

Hutter

The dynamics of avalanches of granular materials from initiation to runout. Part I: analysis. Acta Mech 1991; 86: 201–223.

Hutter

Siegel

Savage

, et al. Two-dimensional spreading of a granular avalanche down an inclined plane part I. Theory. Acta Mech 1993; 100: 37–68.

Bouchut

Westdickenberg

Gravity driven shallow water models for arbitrary topography. Commun Math Sci 2004; 2: 359–389.

Gray

JMNT

Wieland

Hutter

Gravity-driven free surface flow of granular avalanches over complex basal topography. Proc R Soc A Math Phys Eng Sci 1999; 455: 1841–1874.

Pudasaini

Hutter

Rapid shear flows of dry granular masses down curved and twisted channels. J Fluid Mech 2003; 495: 193–208.

Pudasaini

Wang

Hutter

Modelling debris flows down general channels. Nat Hazard Earth Syst Sci 2005; 5: 799–819.

Iverson

Denlinger

RP.

Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory. J Geophys Res Solid Earth 2001; 106: 537–552.

10.

Denlinger

Iverson

RM.

Granular avalanches across irregular three-dimensional terrain: 1. Theory and computation. J Geophys Res Earth Surf 2004; 109: 1–16.

11.

Caleffi

Valiani

Bernini

Fourth-order balanced source term treatment in central WENO schemes for shallow water equations. J Comput Phys 2006; 218: 228–245.

12.

Canestrelli

Siviglia

Dumbser

, et al. Well-balanced high-order centred schemes for non-conservative hyperbolic systems. Applications to shallow water equations with fixed and mobile bed. Adv Water Resour 2009; 32: 834–844.

13.

Castro

Gallardo

Parés

High order finite volume schemes based on reconstruction of states for solving hyperbolic systems with nonconservative products. applications to shallow-water systems. Math Comput 2006; 75: 1103–1135.

14.

Pelanti

Bouchut

Mangeney

, et al. Numerical modeling of two-phase gravitational granular flows with bottom topography. In: Benzoni-Gavage

Serre

(eds) Hyperbolic problems: theory, numerics, applications. Berlin, Heidelberg: Springer, 2008, pp.825–832.

15.

Pelanti

Bouchut

Mangeney

A Roe-type scheme for two-phase shallow granular flows over variable topography. ESAIM Math Model Numer Anal 2008; 42: 851–885.

16.

Pitman

A two-fluid model for avalanche and debris flows. Philos Trans A Math Phys Eng Sci 2005; 363: 1573–1601.

17.

Favreau

Mangeney

Lucas

, et al. Numerical modeling of landquakes. Geophys Res Lett 2010; 37: 1–5.

18.

Kuo

Tai

Bouchut

, et al. Simulation of Tsaoling landslide, Taiwan, based on Saint Venant equations over general topography. Eng Geol 2009; 104: 181–189.

19.

Pirulli

Mangeney

Results of back-analysis of the propagation of rock avalanches as a function of the assumed rheology. Rock Mech Rock Eng 2008; 41: 59–84.

20.

Mangold

Mangeney

Migeon

, et al. Sinuous gullies on Mars: frequency, distribution, and implications for flow properties. J Geophys Res Planets 2010; 115: E11.

21.

Pelanti

Bouchut

Mangeney

A Riemann solver for single-phase and two-phase shallow flow models based on relaxation. Relations with Roe and VFRoe solvers. J Comput Phys 2011; 230: 515–550.

22.

Qamar

Zia

Ashraf

The space-time CE/SE method for solving single and two-phase shallow flow models. Comput Fluids 2014; 96: 136–151.

23.

Zia

Qamar

A kinetic flux-vector splitting method for single-phase and two-phase shallow flows. Comput Math Appl 2014; 67: 1271–1288.

24.

Yan

Ito

Numerical investigation of indoor particulate contaminant transport using the Eulerian-Eulerian and Eulerian-Lagrangian two-phase flow models. Exp Comput Multiphas Flow 2019; 2: 31–40.

25.

Greenberg

Leroux

AY.

A well-balanced scheme for the numerical processing of source terms in hyperbolic equations. SIAM J Numer Anal 1996; 33: 1–16.

26.

Hubbard

Garcia-Navarro

Flux difference splitting and the balancing of source terms and flux gradients. J Comput Phys 2000; 165: 89–125.

27.

Jin

A steady-state capturing method for hyperbolic systems with geometrical source terms. ESAIM Math Model Numer Anal 2001; 35: 631–645.

28.

Kurganov

Petrova

A second-order well-balanced positivity preserving WENO scheme for the Saint-Venant system. Commun Math Sci 2007; 5: 133–160.

29.

Shu

CW.

High-order finite difference and finite volume WENO schemes and discontinuous Galerkin methods for CFD. Int J Comput Fluid Dyn 2003; 17: 107–118.

30.

Harten

Osher

Uniformly high-order accurate nonoscillatory schemes. I. SIAM J Numer Anal 1987; 24: 279–309.

31.

Shu

Osher

Efficient implementation of essentially non-oscillatory shock-capturing schemes. J Comput Phys 1988; 77: 439–471.

32.

Liu

Osher

Chan

Weighted essentially non-oscillatory schemes. J Comput Phys 1994; 115: 200–212.

33.

Jiang

Shu

CW.

Efficient implementation of weighted ENO schemes. J Comput Phys 1996; 126: 202–228.

34.

Shu

CW.

Weighted essentially non-oscillatory schemes on triangular meshes. J Comput Phys 1999; 150: 97–127.

35.

Balsara

Shu

CW.

Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy. J Comput Phys 2000; 160: 405–452.

36.

Xing

Shu

CW.

High-order well-balanced finite difference WENO schemes for a class of hyperbolic systems with source terms. J Sci Comput 2006; 27: 477–494.

37.

Shu

CW.

High order weighted essentially nonoscillatory schemes for convection dominated problems. SIAM Rev 2009; 51: 82–126.

38.

Zhang

Shu

CW.

Positivity-preserving high order finite difference WENO schemes for compressible Euler equations. J Comput Phys 2012; 231: 2245–2258.

39.

Qiu

Hybrid well-balanced WENO schemes with different indicators for shallow water equations. J Sci Comput 2012; 51: 527–559.

40.

Xing

High order finite volume WENO schemes for the Euler equations under gravitational fields. J Comput Phys 2016; 316: 145–163.

41.

Zhu

Qiu

A new type of finite volume WENO schemes for hyperbolic conservation laws. J Sci Comput 2017; 73: 1338–1359.

42.

Zhu

Shu

CW.

A new type of multi-resolution WENO schemes with increasingly higher order of accuracy on triangular meshes. J Comput Phys 2019; 392: 19–33.

43.

Rehman

Ali

Zia

, et al. Well-balanced finite volume multi-resolution schemes for solving the Ripa models. Adv Mech Eng 2021; 13: 1–16.

44.

Ahmed

Rehman

Ali

, et al. A high order multi-resolution WENO numerical scheme for solving viscous quantum hydrodynamic model for semiconductor devices. Results Phys 2021; 23: 104078.

45.

Castro

Macías

Parés

AQ-scheme for a class of systems of coupled conservation laws with source term. Application to a two-layer 1-D shallow water system. ESAIM Math Model Numer Anal 2001; 35: 107–127.

46.

Castro

Garcıa-Rodrıguez

González-Vida

, et al. Numerical simulation of two-layer shallow water flows through channels with irregular geometry. J Comput Phys 2004; 195: 202–235.

47.

Vreugdenhil

CB.

Two-layer shallow-water flow in two dimensions, a numerical study. J Comput Phys 1979; 33: 169–184.

48.

Abgrall

Karni

Two-layer shallow water system: a relaxation approach. SIAM J Sci Comput 2009; 31: 1603–1627.