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

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.


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. 1 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 onedimensional single-phase flow model which was proposed to study avalanches. Later, they extended it to two-dimensions 4 and several authors generalized the results on more complex basal topography. [5][6][7][8][9][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][12][13][14][15] The two-phase shallow granular flow model considered in this article is a variant of Pitman and Le 16 fluid model , comprising of shallow layer of liquid and solid phases. In literature, the impacts of the liquid phase [17][18][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. 20 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][22][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. 20 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 twophase shallow granular flow model. 21 In Yan et al., 24 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, 23 the schemes preserves the positivity of height only for g = 1 2 . Therefore, our motivation was to extend a scheme which can handle different flow conditions for all values of g\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 wellbalanced). 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][26][27][28] In this study, a fifth-order finite volume WENO numerical technique 29 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 30 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 31 developed the finite difference ENO scheme with total variation diminishing Runge-Kutta (TVD-RK) time discritization procedure. Subsequently, Liu et al. 32 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 33 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][35][36][37][38][39][40][41][42][43][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][34][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 r f and r s be the densities of liquid and solid phases of an incompressible flow respectively. The height of the flow be h and f be the solid volume fraction. The other variables are defined as The velocities in x-direction of both solid and fluid phases are denoted by u x s , u x f , respectively. Phase momenta are given by m x s = h x s u x s and m x f = h x f u x f . The continuity and momentum equations for the two-phase shallow system are 21,22 Here, b = def b(x) is basal topography, g is the gravity constant, F D is representing inter-phase drag forces which are given by, where D is a drag function depending generally upon f, ju x f À u x s j and a set of physical parameters s (e.g. particle diameter, specific densities) and g = r f =r s \1. If we take D = 0, then the equations above are somewhat similar to the two-layer shallow flow models. [45][46][47] In this article, we do not consider the drag forces. In compact form the system of equations are expressed as 21 where, The continuity equations for h x s and h x f are conservative, whereas m x s and m x f representing momentum equations for both solid and liquid phases respectively. The system can also be written in following quasi-linear form The conditions on the variables for this model are: , the conditions listed above can be written as follows:, 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 as well as dividing the space O into cells C i = ½x iÀ 1 2 , x i + 1 2 , i = 1, 2, 3:::, N : The center and length of the i À th cell is indicated by in this case by and Dx i respectively. Upon integration of equation (12) over C i , we get where W i (t) = 1 Where, monotone numerical flux is denoted bŷ As a monotone numerical flux, we shall utilize the Lax-Friedrichs flux (LFF), which is defined as follows where q = max W jF 0 (W )j. Cell averages W (x i , t) are approximated by the computational variables W i (t). By using WENO reconstruction, the W À whereŴ l i + 1 2 andW l iÀ 1 2 , for l = 1, 2, 3 are values that have been recreated and described aŝ The nonlinear weights v l andṽ l in equations (16) and (17) and B l = The spatial reconstruction process is now finished. On the right-hand side of equation (13), the nonzero terms is discretized after that, The averaged value at cell C i is gh x s i, while the source fluxes at the cell C i 's interfaces are (b) i 6 1 2 . In addition, for C i 's right interface, we write as The values of cell C i 's left interface and Further, S i denotes the approximation of the expres- Finally, we arrive at the following semi-discrete equation: alternatively, the equation above may be written as Now, we use the third-order TVD RK technique 29 to solve the system of ordinary differential equation (33) as 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 CFL Ã dx max (jl 1 j, jl 2 j, jl 3 j, jl 4 j, jl 5 j) :

Numerical test problems
To validate the proposed WENO scheme several numerical test problems are considered here. In all the experiments, we set g\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. 22 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.
This problem is also discussed in References. 15,22,23 The interface is located at x = 0 initially. The values of h, f, u x s , u x f at right and left of interface are given below (h, f, u x s , u x f ) R = (2, 0:4, À 0:9, À 0:1): Here, we take b(x) = constant, g = 9:81 and L = 10.
The domain À5 ł x ł 5 is discritized into 300 grid points. At t = 0:5 the results for h, h x s , h x f , the solid volume fraction f and the phase velocities u x s , u x f are shown in Figure 1. The approximated solution consists of 1 -rarefaction wave, 2 -shock wave, and 3rarefaction 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).
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: 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 f = 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 À f) vanishes. In addition, four more rarefaction waves and one shock wave are observed in a pure solid area where (h x f = 0) is observed. One can also observe that results are free of oscillations.
Here, b(x) = constant and g = 1:0. In this test case we take h = 0 below the tolerance 2 = 10 À5 . The domain À10 ł x ł 10 is divided into 300 cells. Though inter-phase drag forces are not considered here, still ju x s À u x f j ! 0 as h ! 0. The results are depicted in the Figure 3. The hyperbolic conditions are maintained as h ! 0. Also a noticeable change in of f (solid volume fraction) can also be observed. 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 Here, g = 9:81 and g = 0:98. The initial interface is located at x = 0:5. The domain 0 ł x ł 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.
This test problem is related to the formation of a dry bed zone. The left and right state data for this problem is, (h, f, u x s , u x f ) R = (0:1, 0:7, 3:0, 3:0), if x.0 : The initial discontinuity is at x = 0 and b(x) = constant. The gravitational constant g = 9:81 and g = 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.

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 g\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.

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.