High order multi-resolution WENO scheme with AUSMV numerical flux for solving the two-phase flows

This article presents the development of a fifth-order multi-resolution finite volume weighted essentially non-oscillatory (WENO) scheme combined with the advection upstream splitting method based on flux vector splitting (AUSMV) numerical flux for analyzing two-phase flow in both horizontal and vertical pipelines. The drift flux flow model comprises of two separate mass conservation equations for each phase for liquid and gas and one momentum equation for mixture and submodels for thermodynamics and hydrodynamics. The two mass conservation equations describe the behavior of each phase in the flow. The mixture-momentum equation takes into account the frictional and gravitational forces acting on the mixture of both phases. The thermodynamic and hydrodynamic submodels provide additional information to fully describe the flow and close the drift flux model. In the presence of these source terms and submodels, it is a challenging task to develop a high order efficient and accurate numerical schemes. The proposed numerical technique captures the peaks of pressure wave, suppresses the erroneous oscillations at the transition zones and resolves the discontinuities more efficiently and accurately. The accuracy of proposed numerical technique is verified by solving the various test problems. Furthermore, the solution obtained by developed numerical technique are compared to those attained with the high-resolution improved CUP and simple finite volume WENO numerical schemes.


Introduction
Two-phase flow models are important because they help in predicting the behavior of fluids in various industrial processes, such as oil and gas production, power generation, and cooling systems.2][3] Two-phase flow models play a critical role in ensuring the safe, efficient, and cost-effective operation of many industrial processes.To study the behavior of two-phase flows, researchers have developed several two-phase flow models and a wide range of numerical approaches.Despite their importance, these models can be difficult to analyze numerically due to the complexities of the interactions between the two phases.More precisely, the two-phase flow models offer various numerical difficulties, such as the discontinous solutions, formation of shocks, and loss of hyperbolicity.To overcome such type of numerical difficulties, Zuber and Findly 4 introduced the drift flux model.Later, various authors improved and modified the drift flux flow model, for example see the articles [5][6][7][8][9] and references therein.In the drift flux flow model, the one momentum equation for mixture is utilized instead of using separate equations of momentum and introduced the hydrodynamic submodels.This approach helps in eliminating certain problematic terms associated with phase interactions.The use of a mixture momentum equation provides a more accurate representation of the fluid flow.Further, the hydrodynamic submodels help in accurately computing the phase velocities in the solution domain.This model provides an efficient means of understanding multi-phase flows in well-bores.
There have been wide range of studies in the literature in which researchers have developed numerical schemes to approximate and analyze the drift flux model in various flow systems such as pipelines, well-bores, and marine risers.Initially, the authors developed the Roe-type numerical scheme see for example, Romate. 5he two-phase flow model is then further analyzed using AUSM and relaxation type numerical techniques, highresolution hybrid-upwind, hybrid flux splitting, for detail see Fjelde and Karlsen, 6 Steinar and Fjelde, 7 Baudin et al., 8 and Evje and Fjelde. 10,11Futhermore, weakly implicit-mixture flux (WIMF), multi stage approach (MUSTA) based centered schemes, space-time CESE schemes 9,12,13 were developed for finding the solutions of model considered in this article.Newtonian and non-Newtonian experimental methods are used to explore both horizontal and inclined pipelines by using the drift flux flow model. 14Next, the two-phase flow in horizontal and vertical channels are simulated in the article. 15All above mentioned numerical schemes are of first or second order and mostly the research articles discussed the two-phase flow only for horizontal pipe-lines.Recently, the third order and fifth order WENO schemes with the central upwind flux 16 are developed for the drift flux flow model.In this article, 16 the authors considered the drift flux model without source terms and only tested the proposed numerical schemes for two-phase flow in horizontal pipe lines.
The main purpose of this study is to develop the fifth order finite volume multi-resolution WENO scheme 17 with AUSMV flux 18 for analyzing the behavior of twophase flow in both vertical and horizontal pipelines, 10 with and with out source terms.This new hybrid numerical scheme is more efficient, accurate, and less dissipative than the simple WENO and improved central upwind numerical schemes.The main idea of this hybrid scheme is to calculate the numerical flux by using the AUSMV approach which is a modified version of the original AUSM numerical flux that was developed to handle both inviscid and viscous flows.It uses the concept of splitting the total flux into convective flux and pressure flux and includes an artificial viscosity term to damp the numerical oscillations.Further, the cell interface values are approximated by multi-resolution WENO reconstruction which is the extension of WENO numerical scheme.The multiresolution WENO method uses a combination of coarse and fine grids to represent the solution and employs a multi-level correction process to improve the accuracy of the solution.The method uses a local adaptive grid refinement strategy, which allows for a more accurate representation of the solution in regions where high gradients or discontinuities are present.
1][22] In this type of the finite volume WENO scheme, they used the convex combination of all candidate stencils and increased the accuracy order.After that, these numerical methods have been developed, updated, and expanded for several disciplines of engineering and science for further information, see Hu and Shu, 23 Balsara and Shu, 24 Shu, 25 Xing and Shu, 26 Shu, 27 Zhang and Shu, 28 and Zhu and Qiu. 29Later, the authors 17,30 developed multi-resolution finite volume and finite difference WENO schemes for solving hyperbolic conservation laws.On the other hand, Initially, the AUSM numerical scheme was introduced by Liou and Steffen for the Euler equation. 31Wada and Liou later proposed an improved version of the scheme known as AUSMDV in Wada and Liou. 32The authors of the article 18 then introduced a more precise and reliable numerical scheme called AUSMV, which combines the flux vector-splitting (FVS) and AUSM methods.For more information, please see the references Evje and Fjelde, 18 Liou and Steffen, 31 and Wada and Liou 32 and the sources cited therein.
The rest of the article is organized as follows.In Section ''Drift flux model for two-phase flows,'' a twophase drift flux flow model is described.In section ''Construction of multi-resolution WENO scheme along with AUSMV numerical flux for the drift flux model,'' drift flux flow model is derived by the hybrid numerical scheme (multi-resolution WENO along with AUSMV flux).Test problems are given in section ''Test problems.''Also, the solution profiles are compared with those attained from improved central-upwind and simple WENO techniques 33 respectively.The conclusions are presented in last section.

Drift flux model for two-phase flows
The drift flux flow model for two phase-flows is mathematically expressed and eigenvalues are explained in this part.As discussed earlier, the drift flux flow model is composed of three partial differential equations, two seprate mass conservation equations for both phases for liquid and gas and one mixture momentum equation.These equations are mathematically expressed as follows where u l , u g denote the liquid and gas velocity, a l , a g represents the volume fraction of liquid and gas r l , and r g for liquid and gas density, respectively.The mutual pressure is represented by p and the source term is denoted by s.
where, g = 9:81 m=s 2 and the inner diameter of pipeline is denoted by D. The average velocity term u mix and viscosity term m mix in the frictional forces F f are Further, three differential equations of drift flux model (1) contain seven unknown variables a l , a g ,r l , r g , u l , u g , and p, therefore, for obtaining the solution, four additional relations are required, which are and where a l , a g , r l, 0 , p l, 0 , C o , and u d represent the velocities of sound in liquid and gas phases, reference liquid density, reference liquid pressure, flow dependent parameter, and drift gas velocity.The values of these parameter given in the Table 1.

Eigenvalues calculation
The system of equations ( 1) are expressed in the compact form, given as where l +a g r g u 2 g +p) and W=(w 1 , w 2 , w 3 ) =(a l r l , a g r g , a l r l u l +a g r g u g ).Now, the equation ( 8) can be expressed as It is observed that the system indicates the flux vector cannot be fully expressed in terms of the conserved variables.As a result, the exact value of the sound velocity cannot be determined.In the articles, 6,7 the authors have discussed the full detail for approximating the value of sound velocity.we have omited the expanation of procedure here, but only describe the eigenvalues.Hence, the eigenvalues for the given model is Construction of multi-resolution WENO scheme along with AUSMV numerical flux for the drift flux model In this part, two-phase drift flux flow model is derived by fifth order multi-resolution WENO scheme along with AUSMV numerical flux.This model are expressed in the compact form, given as Here F(W) is the total flux, which is the sum of convective flux F c and pressure flux F p , that is expressed as follows further, the convective flux is also the sum of two fluxes, that is F c, l (W) is the liquid convective flux and F c, g (W) is the gas convective flux expressed as For the construction of hybrid numerical scheme, firstly, the given domain is discretized into cells, given as C i = x iÀ 1 2 , x i + 1 2 h i with, i = 1, :::, N : The length and center of iÀth cell is defined as Dx i : = x iÀ 1 2 + x i + 1 2 and respectively.Further, integrating the equation (12) Conservative scheme is used to approximate the equation ( 15), described as follows Here, the numerical flux is reprsented by 6 are pointwise approximations to W x i + 1 2 , t .Generally, the numerical flux is expressed as where, F g, L = (0, 1, u g, L ) T and F l, L = (1, 0, u l, L ) T .
Here, the interface mass fluxes that can be selected in different ways are (a l r l u l ) i + 1 2 , and (a g r g u g ) i + 1 2 , so we are using the AUSMV numerical flux, that is expressed as ) and the viscous term D q, i + 1 2 is defined as ju q, i + 1 2 j ((a q r q ) R À (a q r q ) L ).Further, we defined the velocity splitting function denoted as Ṽ 6 is where, the expressions for X L and X R are given as follows Next, we discretized the pressure term p i + 1 2 as follows further, we defined the pressure splitting function, Finally, the speed of sound in two-phase mixture c i + 1 2 , is explained at the cell-interface here, and E is defined as c(a g ) = a l , for a g \E, a, for E ł a g ł 1 À E a g , for a g .E: Parameter E is taken as 0:001 for smooth transition of sound velocity from two-phase to single-phase.Further, the point-wise approximations and are determined through the linked cell average values W i by WENO reconstruction.Three central candidate stencils are chosen for the fifth order WENO reconstruction.These stenciles are further defined as, g, and then reconstruct the zeroth, second, and fourth degree polynomials p 1 (x),p 2 (x), and p 3 (x) respectively, that satisfies x jÀ 1 2 p 3 (x)dx=W j , here, j=i À 2,i À 1,i,i+1,i+2: To be more specific, these polynomials are explicitly defined as p 1 (x), p 2 (x), and p 3 (x) The reconstructed values point-wisely are attained by the following relations where , are reconstructed values for l = 1, 2, 3, and expressed as with Where j 1, 2 + j 2, 2 = 1, j 1, 3 + j 2, 3 + j 3, 3 = 1, and j 2, 2 , j 3;3 6 ¼ 0. The linear weights are denoted by j 1, 2 , j 2, 2 in these formulations.We fix the linear weights as to measure the accuracy in smooth regions and in non-smooth regions maintain a balancing between the sharp and essentially non-oscillatory shock transitions, j 1, 2 = 1=11, j 2, 2 = 10=11, j 1, 3 = 1=111, j 2, 3 = 10=111, and j 3, 3 = 100=111 for the fifth order approximation, that is explained in Zhu and Shu. 17 While, non linear weights are v l given in equation ( 28) are expressed as follows Here, j l represent the linear weights and B l represent smoothness indicators and in all the simulations e is taken as 10 À10 .While expression for given smoothness indicators in general form is expressed as follows Here, for m = 2, 3, k is found by the formula k = 2(m À 1).The values of B 1 , B 2 , and B 3 are defined in the same way as given in Zhu and Shu. 17 To be more accurate, these smooth indicators are expressed as follow where To resolve the discontinuities efficiently, B 1 is expressed as follows In equation ( 47), To avoid the denominator to become zero e is taken as a small positive number.Finally, we obtained The term t in equation ( 36) is stated as This completes the derivation of WENO technique with AUSMV numerical flux for the two-phase drift flux flow model.Next, we will approximate the source term in equation (8).The approximation of the term x iÀ 1 2 S(W)dx is denoted by S i and given as follows where t denotes the transpose.Thus, we left the system of ODEs which are the above system can also be rewritten as Now, the third order TVD RK method 17 is used to compute the system equation ( 52), as follow (W (1) + dtL(W (1) )), where, by using the relation CFL Ã dx max (jl 1 j, jl 2 j, jl 3 j) time step dt is determined, where CFL denotes the Courant-Friedrichs-Lewy condition.L(W) is the spatial operator in the above equation.

Test problems
In this part, we have compared results of high resolution numerical schemes such as improved central upwind 33 with our developed hybrid scheme (MR WENO).Further, we have compared the results of developed scheme and WENO scheme with Lax-Friedrichs numerical flux (WENO) for the second test problem.
Test Problem 1: This test problems is selected from Kuila et al. 34 for measuring the efficiency of developed technique.Here, no-slip condition is used, it means u = u g = u l and also set the source term s = 0.The initial data is given as (r g , r l , u, a g , a l )(x, 0) = (1, 1, 0, 500, 100), x.0, (6, 1, 0, 500, 100), & For test problem 1, p, a the pressure and speed of sound are determined by using the following relations ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p r g a g + r l a l r and a 2 g r g + a 2 l r l : The length of the pipe is considered 800 m long.In Figure 1 at t = 0:5 s and on 400 grid cells we have calculated the numerical results for the gas and liquid density and also computed results for pressure and velocity.The solution domain consists of shock wave, rarefaction wave, and contact wave.The developed and improved central upwind techniques have good agreement with each other but we observed that the less diffusivity occurs in our developed numerical scheme rather then the above mentioned scheme.
Test problem 2: To verify the accuracy and effectiveness of the suggested technique we have selected this problem from Evje and Fjelde. 10Moreover, a length of 100 m horizontal channel is taken into consideration together with the gas-slip relation (7) in algebraic form.Flow dependent parameter C o = 1:07 is taken and drift velocity in the gas-slip relation are taken as u d = 0:216 m=s.Now, Initial data for this problem is expressed as, (a g , u l , p)(x, 0) = (0:55, 10:37, 80450), x\50, (0:55, 0:561, 24282), x.50:

&
Hybrid numerical scheme (MR WENO) are used to determined the profiled solution for gas velocity, liquid velocity, and pressure volume fraction, simple WENO, and improved central upwind numerical scheme (MCUP) at time t = 1:0 s, shown in the Figure 2. The reference solution is computed on 8000 grid points.Domain of the solution consists of two shock waves and one rare fraction waves.By using reference solution, the L 1 À error is computed, given in Table 2.This table and solution profiles show that the MR WENO numerical scheme is more efficient then the fifth order WENO and improved CUP numerical schemes.are resolved by derived numerical scheme (MR WENO) more accurately and efficiently then the highresolution numerical scheme (MCUP).
Test problem 4: we have taken this problem from the article. 6This test problem is more challenging for the numerical schemes since this test problem, considered frictional forces (2).Initially, in the horizontal pipe we have considered 10 5 Pa pressure whose length is 1000 m having diameter 0:1 m.In start dormant liquid is filled in the pipe 750 m and in remaining portion 99 % liquid and only 1 % of gas.Further, for 0:0025 seconds, the liquid at the rate 0-0:3 kg=s is injected to the inlet of the pipe.The profiled solution of liquid velocities at times t = 0:4s and t = 1:3s are shown in top row of the Figure 4. While, the profiled solution for the pressure at times t = 0:4s and t = 1:3s are shown in the bottom row of the Figure 4. Here, the further detail about the physics of solution profiles clearly shows that the comparison of improved central upwind numerical scheme and the WENO numerical scheme of liquid velocity profiled solution for a time t = 1:3 s clearfies that WENO is much better then mentioned schemes.Test problem 5: This problem is deemed to be the most difficult challenge for numerical methods taken from Fjelde and Karlsen. 6The source term in    Afterward, to increase the mass flow rates of the gas and liquid from 0:02 kg=s to 3:0 kg=s respectively,over the course of the first 10 s gas and liquid are introduced into the pipe inlet at L = 0 m.The pressure is kept constant 10 5 pa at the pipe outlet at L = 1000 m.In the Figure 5, the top row of figures depict the liquid profile mass flow rates at different positions x = 250m, x = 750m, respectively, and the bottom row of figures show gas mass flow rates profiles at different position x = 250m, x = 750m within the solution space.Our final results are in close agreement with those presented in Fjelde and Karlsen, 6 therefore, a detailed explanation of these profiles is omitted.For further information on the physical meaning behind these profiles, the reader is directed to Fjelde and Karlsen. 6The profiled solution, obtained by derive hybrid numerical scheme (MR WENO) and improved central upwind scheme (MCUP), in the Figure 5 clearly depict that rather then MCUP numerical scheme MR WENO scheme has better ability to capture the sharp peaks.
Test problem 6: Test problem 6 is taken to analyze the two-phase flow in the vertical marine riser of length 3048 m.In off-shore drilling, the riser yields the communication and circulation capability between the well-head in the sea floor and surface. 35n this problem, to maintain the hydrostatic pressure in the drilling process water based mud is used, for more detail the reader is referred to Velmurugan. 35h Hence, the water and methane are used for the liquid and gas phases respectively.In this case, both terms of the source terms are considered and the drift velocity in the gas slip relation 7 is defined as follows This type of drift velocity 35 used for bubble flow in the slug flow.Further, the choke pressure is considered 100 bars.This pressure enforces the pressure gradient in the system, more precisely this pressure helps to prevent dissipation of pressure in the system.The vertical marine riser is discretized into 100 grid points.The solution profiles shown in the Figure 6 are obtained at t = 250 s.The pressure decreases p is mainly due to the frictional or viscous F w and gravitational F g losses, the decreasing pressure causes the expansion in gas. Threfore the gas volume fraction and gas phase velocities increases as shown in the Figure 6.

Conclusions
In this research, a high order hybrid numerical scheme is developed for simulating the two-phase flow in horizontal and vertical channels, more precisely, a fifthorder multi-resolution WENO is utilized for the spatial reconstruction, and the AUSMV numerical scheme is used to the calculate the fluxes of conserved quantities across the computational cells in a discretized domain.This new hybrid numerical scheme effectively resolves contact discontinuities, conserves the positivity of flow variables, and eliminates spurious oscillations in transition zones or in the vicinity of sharp gradients.A comprehensive set of test cases demonstrates the precision and dependability of the newly developed technique.Additionally, the profiled solution attained by new method are compared to those produced by the improved CUP and WENO with Lax-Friedrichs flux (simple WENO) numerical approaches.The results show that our developed numerical technique is robust, more efficient, and less diffusivity as compared to the improved central upwind and simple WENO numerical methods.

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.

Test problem 3 :
Domain of the solution consits of two rare-faction waves and one contact wave in this problem.The length of the horizontal pipe, final time, and values of the parameters involved in the gas-slip relation are considered same as given in the previous test problem.The data of the problem is expressed as follows (a g , u l , p)(x, 0) = (0:35, 1:868, 192170), x\50, (0:30, 14:47, 196690), x.50: & The profiled solution obtained by MR WENO and MCUP numerical schemes are shown in the Figure 3. Precisely, the profiled solution consist of contact waves and rare-faction.Moreover, the contact discontinuities

Figure 3 .
Figure 3.The solution profiles of test problem 3 at t = 1:0 s.

Figure 4 . 5
Figure 4.The solution profiles of pressure pulse and liquid velocity at different times.

Figure 5 .
Figure 5.The solution profiles of test problem 5.