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

Abstract

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.

Keywords

Multi-resolution two-phase flow AUSMV flux WENO scheme drift flux flow model source terms shock wave

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. Further, these models are extensively used for: understanding fluid behavior in pipelines, heat exchangers, and other equipment; designing and optimizing industrial processes; reducing operational costs and improving efficiency and preventing safety hazards, such as explosions and leaks.^1–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⁴ introduced the drift flux model. Later, various authors improved and modified the drift flux flow model, for example see the articles^5–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.⁵ The two-phase flow model is then further analyzed using AUSM and relaxation type numerical techniques, high-resolution hybrid-upwind, hybrid flux splitting, for detail see Fjelde and Karlsen,⁶ Steinar and Fjelde,⁷ Baudin et al.,⁸ and Evje and Fjelde.^10,11 Futhermore, 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.¹⁴ Next, the two-phase flow in horizontal and vertical channels are simulated in the article.¹⁵ All 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¹⁶ are developed for the drift flux flow model. In this article,¹⁶ 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¹⁷ with AUSMV flux¹⁸ for analyzing the behavior of two-phase flow in both vertical and horizontal pipelines,¹⁰ 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 multi-resolution 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.

Initially, in 1994, Liu et al.¹⁹ implemented the finite volume WENO numerical scheme of third order, which was the refined form of the essentially non-oscillatory numerical approach.^20–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,²³ Balsara and Shu,²⁴ Shu,²⁵ Xing and Shu,²⁶ Shu,²⁷ Zhang and Shu,²⁸ and Zhu and Qiu.²⁹ Later, 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.³¹ Wada and Liou later proposed an improved version of the scheme known as AUSMDV in Wada and Liou.³² The authors of the article¹⁸ 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,¹⁸ Liou and Steffen,³¹ and Wada and Liou³² and the sources cited therein.

The rest of the article is organized as follows. In Section “Drift flux model for two-phase flows,” a two-phase 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³³ 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

\begin{matrix} \frac{\partial}{\partial_{t}} (a_{l} ρ_{l}) + \frac{\partial}{\partial_{x}} (a_{l} ρ_{l} u_{l}) = 0, \\ \frac{\partial}{\partial_{t}} (a_{g} ρ_{g}) + \frac{\partial}{\partial_{x}} (a_{g} ρ_{g} u_{g}) = 0, \\ \frac{\partial}{\partial_{t}} (a_{l} ρ_{l} u_{l} + a_{g} ρ_{g} u_{g}) + \frac{\partial}{\partial_{x}} (a_{l} ρ_{l} u_{l}^{2} + a_{g} ρ_{g} u_{g}^{2} + p) = s, \end{matrix}

(1)

where $u_{l}$ , $u_{g}$ denote the liquid and gas velocity, $a_{l}$ , $a_{g}$ represents the volume fraction of liquid and gas $ρ_{l}$ , and $ρ_{g}$ for liquid and gas density, respectively. The mutual pressure is represented by $p$ and the source term is denoted by $s$ .

s = F_{f} + G_{f} = - \frac{32 u_{mix} μ_{mix}}{D^{2}} + (a_{g} ρ_{g} + a_{l} ρ_{l}) g,

(2)

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 $μ_{mix}$ in the frictional forces $F_{f}$ are

μ_{mix} = a_{l} μ_{l} + a_{g} μ_{g} and u_{mix} = a_{g} u_{g} + a_{l} u_{l} .

(3)

Further, three differential equations of drift flux model (1) contain seven unknown variables $a_{l}, a_{g}$ , $ρ_{l}, ρ_{g}$ , $u_{l}, u_{g}$ , and $p$ , therefore, for obtaining the solution, four additional relations are required, which are

ρ_{l} = ρ_{l, 0} + \frac{p - p_{l, 0}}{α_{l}^{2}},

(4)

ρ_{g} = \frac{p}{α_{g}^{2}},

(5)

a_{l} + a_{g} = 1 .

(6)

and

u_{g} = C_{o} u_{mix} + u_{d},

(7)

where $α_{l}$ , $α_{g}$ , $ρ_{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.

Table 1.

Parameter values.

Parameters	Value	Description
$μ_{g}$	$5 . 10^{- 6} Pa . s$	Gas phase viscosity
$μ_{l}$	$5 . 10^{- 2} Pa . s$	Liquid phase viscosity
$α_{g}$	$316 m / second$	Velocity of sound in gas phase (water)
$α_{l}$	$1000 m / second$	Velocity of sound in liquid phase
$α_{g}$	$446 m / second$	Velocity of sound in methane gas phase (for Test problem 6)
$ρ_{l, 0}$	$1000 Kg / m^{3}$	Reference density liquid
$p_{l, 0}$	$10^{5} Pa$	Reference pressure liquid
$p_{choke}$	$100 bars$	Choke pressure

Eigenvalues calculation

The system of equations (1) are expressed in the compact form, given as

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

(8)

where

$W = (\begin{matrix} a_{l} ρ_{l} \\ a_{g} ρ_{g} \\ a_{l} ρ_{l} u_{l} + a_{g} ρ_{g} u_{g} \end{matrix}), F (W) = (\begin{matrix} a_{l} ρ_{l} u_{l} \\ a_{g} ρ_{g} u_{g} \\ a_{l} ρ_{l} u_{l}^{2} + a_{g} ρ_{g} u_{g}^{2} + p \end{matrix}), S (W) = (\begin{matrix} 0 \\ 0 \\ s \end{matrix})$ ,

with, $S (W) = (s_{1}, s_{2}, s_{3}) = (0, 0, s)$ , $F (W) = (f_{1}, f_{2}, f_{3}) = (a_{l} ρ_{l} u_{l}, a_{g} ρ_{g} u_{g}, a_{l} ρ_{l} u_{l}^{2} + a_{g} ρ_{g} u_{g}^{2} + p)$ and $W = (w_{1}, w_{2}, w_{3}) = (a_{l} ρ_{l}, a_{g} ρ_{g}, a_{l} ρ_{l} u_{l} + a_{g} ρ_{g} u_{g})$ . Now, the equation (8) can be expressed as

\partial_{t} (\begin{matrix} w_{1} \\ w_{2} \\ w_{3} \end{matrix}) + \partial_{x} (\begin{matrix} w_{1} w_{l} \\ w_{2} w_{g} \\ w_{1} w_{l}^{2} + w_{2} w_{g}^{2} + p (w_{1}, w_{2}) \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ s \end{matrix}) .

(9)

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

λ_{1} = u - α,

λ_{2} = u,

(10)

λ_{3} = u + α .

With

α = \sqrt{\frac{p}{a_{g} ρ_{l} (1 - C_{o} a_{g})}} .

(11)

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

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

(12)

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

F (W) = F_{c} (W) + F_{p} (W) = (\begin{matrix} a_{l} ρ_{l} u_{l} \\ a_{g} ρ_{g} u_{g} \\ a_{l} ρ_{l} u_{l}^{2} + a_{g} ρ_{g} u_{g}^{2} \end{matrix}) + (\begin{matrix} 0 \\ 0 \\ p \end{matrix}),

(13)

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

\begin{matrix} F_{c} (W) = F_{c, l} (W) + F_{c, g} (W) & = (\begin{matrix} a_{l} ρ_{l} u_{l} \\ 0 \\ a_{l} ρ_{l} u_{l}^{2} \end{matrix}) + (\begin{matrix} 0 \\ a_{g} ρ_{g} u_{g} \\ a_{g} ρ_{g} u_{g}^{2} \end{matrix}) \\ = a_{l} ρ_{l} u_{l} (\begin{matrix} 1 \\ 0 \\ u_{l} \end{matrix}) + a_{g} ρ_{g} u_{g} (\begin{matrix} 0 \\ 1 \\ u_{g} \end{matrix}) \end{matrix}

(14)

For the construction of hybrid numerical scheme, firstly, the given domain is discretized into cells, given as $C_{i} = [x_{i - \frac{1}{2}}, x_{i + \frac{1}{2}}]$ with, $i = 1, . . ., N .$ The length and center of $i - th$ cell is defined as $Δ x_{i} . = x_{i - \frac{1}{2}} + x_{i + \frac{1}{2}}$ and $x_{i} = \frac{(x_{i - \frac{1}{2}} + x_{i + \frac{1}{2}})}{2}$ respectively. Further, integrating the equation (12) over the $i - th$ cell $C_{i}$

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

(15)

where, $\bar{W} (x_{i}, t) = \frac{1}{Δ x_{i}} \int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} W (x, t) dx$ . Conservative scheme is used to approximate the equation (15), described as follows

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

(16)

Here, the numerical flux is reprsented by ${\hat{F}}_{i + \frac{1}{2}} = F_{i + \frac{1}{2}} (W^{-}, W^{+})$ and, $W_{i + \frac{1}{2}}^{\pm}$ are pointwise approximations to $W (x_{i + \frac{1}{2}}, t)$ . Generally, the numerical flux is expressed as

\begin{matrix} F_{i + \frac{1}{2}} (W^{-}, W^{+}) = \\ \frac{1}{2} ({(a_{l} ρ_{l} u_{l})}_{i + \frac{1}{2}} (F_{l, L} + F_{l, R}) - | a_{l} ρ_{l} u_{l})_{i + \frac{1}{2}} | (F_{l, R} - F_{l, L})) + \frac{1}{2} \\ ({(a_{g} ρ_{g} u_{g})}_{i + \frac{1}{2}} (F_{g, L} + F_{g, R}) - | a_{g} ρ_{g} u_{g})_{i + \frac{1}{2}} | (F_{g, R} - F_{g, L})) + p_{i + \frac{1}{2}}, \end{matrix}

(17)

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} ρ_{l} u_{l})_{i + \frac{1}{2}}$ , and $(a_{g} ρ_{g} u_{g})_{i + \frac{1}{2}}$ , so we are using the AUSMV numerical flux, that is expressed as

\begin{matrix} (a_{q} ρ_{q} u_{q})_{i + \frac{1}{2}} = u_{q, i + \frac{1}{2}} \frac{1}{2} ((a_{q} ρ_{q})_{L} + (a_{q} ρ_{q})_{R}) \\ - \frac{1}{2} D_{q, i + \frac{1}{2}}, where q = l, g, \end{matrix}

(18)

where $u_{q, i + \frac{1}{2}} = {\tilde{V}}^{+} (u_{q, L}, c_{i + \frac{1}{2}}, X_{L}) + {\tilde{V}}^{-} (u_{q, R}, c_{i + \frac{1}{2}}, X_{R})$ and the viscous term $D_{q, i + \frac{1}{2}}$ is defined as $| u_{q, i + \frac{1}{2}} | ((a_{q} ρ_{q})_{R} - (a_{q} ρ_{q})_{L})$ . Further, we defined the velocity splitting function denoted as ${\tilde{V}}^{\pm}$ is

\begin{matrix} {\tilde{V}}^{\pm} (u, c, X) = \\ {\begin{matrix} X V^{\pm} (u, c) + (1 - X) \frac{1}{2} (u \pm | u |), & | u | < c, \\ \frac{1}{2} (u \pm | u |), & otherwise . \end{matrix} \end{matrix}

(19)

Next, $V$ is expressed as follows

\begin{matrix} V^{\pm} (u, c) = {\begin{matrix} \pm \frac{1}{4} {(u \pm c)}^{2}, & | u | < c, \\ \frac{1}{2} (u \pm | u |), & otherwise, \end{matrix} \end{matrix}

(20)

where, the expressions for $X_{L}$ and $X_{R}$ are given as follows

X_{L} = \frac{2 {(\frac{ρ_{q}}{α_{q}})}_{L}}{{(\frac{ρ_{q}}{α_{q}})}_{L} + {(\frac{ρ_{q}}{α_{q}})}_{R}}, X_{R} = \frac{2 {(\frac{ρ_{q}}{α_{q}})}_{R}}{{(\frac{ρ_{q}}{α_{q}})}_{L} + {(\frac{ρ_{q}}{α_{q}})}_{R}} .

(21)

Next, we discretized the pressure term $p_{i + \frac{1}{2}}$ as follows

p_{i + \frac{1}{2}} = P^{+} (u_{mix, L}, c_{i + \frac{1}{2}}) p_{L} + P^{-} (u_{mix, R}, c_{i + \frac{1}{2}}) p_{R},

(22)

further, we defined the pressure splitting function,

P^{\pm} (u_{mix}, c) = V^{\pm} (u_{mix}, c) {\begin{matrix} \frac{1}{c} (\pm 2 - \frac{u_{mix}}{c}), & | u_{mix} | < c, \\ \frac{1}{u_{mix}}, & otherwise . \end{matrix}

(23)

Finally, the speed of sound in two-phase mixture $c_{i + \frac{1}{2}}$ , is explained at the cell-interface here,

c_{i + \frac{1}{2}} = max (c_{L}, c_{R}),

and $ϵ$ is defined as

c (a_{g}) = {\begin{matrix} α_{l}, & for a_{g} < ϵ, \\ α, & for ϵ \leq a_{g} \leq 1 - ϵ \\ α_{g}, & for a_{g} > ϵ . \end{matrix}

(24)

Parameter $ϵ$ is taken as $0.001$ for smooth transition of sound velocity from two-phase to single-phase. Further, the point-wise approximations $W_{i + \frac{1}{2}}^{-}$ and $W_{i + \frac{1}{2}}^{+}$ are determined through the linked cell average values ${\bar{W}}_{i}$ by WENO reconstruction. Three central candidate stencils are chosen for the fifth order WENO reconstruction. These stenciles are further defined as, $S_{1} (i) = {C_{i}}, S_{2} (i) = {C_{i - 1}, C_{i}, C_{i + 1}}$ , $S_{3} (i) = {C_{i - 2}, C_{i - 1}, C_{i}, C_{i + 1}, C_{i + 2}}$ , and then reconstruct the zeroth, second, and fourth degree polynomials $p_{1} (x), p_{2} (x)$ , and $p_{3} (x)$ respectively, that satisfies

\frac{1}{Δ x_{i}} \int_{x_{j - \frac{1}{2}}}^{x_{j + \frac{1}{2}}} p_{1} (x) dx = {\bar{W}}_{j}, here, j = i,

(25)

\frac{1}{Δ x_{i}} \int_{x_{j - \frac{1}{2}}}^{x_{j + \frac{1}{2}}} p_{2} (x) dx = {\bar{W}}_{j}, here, j = i - 1, i, i + 1,

(26)

\frac{1}{Δ x_{i}} \int_{x_{j - \frac{1}{2}}}^{x_{j + \frac{1}{2}}} p_{3} (x) dx = {\bar{W}}_{j}, here, j = i - 2, i - 1, i, i + 1, i + 2 .

(27)

To be more specific, these polynomials are explicitly defined as $p_{1} (x), p_{2} (x)$ , and $p_{3} (x)$

\begin{matrix} p_{1} (x) = {\bar{W}}_{i}, \\ p_{2} (x) = \frac{{\bar{W}}_{i + 1} - 2 {\bar{W}}_{i} + {\bar{W}}_{i - 1}}{2 (Δ x_{i})^{2}} {(x - x_{i})}^{2} + \frac{{\bar{W}}_{i + 1} - {\bar{W}}_{i - 1}}{2 (Δ x_{i})} (x - x_{i}) \\ + \frac{- {\bar{W}}_{i + 1} + 26 {\bar{W}}_{i} - {\bar{W}}_{i - 1}}{24}, \\ p_{3} (x) = \frac{1}{1920} [\frac{(80 {\bar{W}}_{i - 2} - 320 {\bar{W}}_{i - 1} + 480 {\bar{W}}_{i} - 320 {\bar{W}}_{i + 1} + 80 {\bar{W}}_{i + 2})}{{(Δ x_{i})}^{4}} {(x - x_{i})}^{4} \\ + \frac{(160 {\bar{W}}_{i - 2} - 320 {\bar{W}}_{i - 1} + 320 {\bar{W}}_{i + 1} - 160 {\bar{W}}_{i + 2})}{{(Δ x_{i})}^{3}} {(x - x_{i})}^{3} \\ + \frac{(120 {\bar{W}}_{i - 2} - 1440 {\bar{W}}_{i - 1} + 2640 {\bar{W}}_{i} - 1440 {\bar{W}}_{i + 1} + 120 {\bar{W}}_{i + 2})}{{(Δ x_{i})}^{2}} {(x - x_{i})}^{2} \\ + \frac{(200 {\bar{W}}_{i - 2} - 1360 {\bar{W}}_{i - 1} + 1360 {\bar{W}}_{i + 1} - 200 {\bar{W}}_{i + 2})}{(Δ x_{i})} (x - x_{i}) \\ + (9 {\bar{W}}_{i - 2} - 116 {\bar{W}}_{i - 1} + 2134 {\bar{W}}_{i} - 116 {\bar{W}}_{i + 1} + 9 {\bar{W}}_{i + 2})] . \end{matrix}

The reconstructed values point-wisely $W_{i + \frac{1}{2}}^{+}$ and $W_{i + \frac{1}{2}}^{-}$ are attained by the following relations

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},

(28)

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},

(29)

where ${\hat{W}}_{i + \frac{1}{2}}^{l}$ and ${\tilde{W}}_{i - \frac{1}{2}}^{l}$ , are reconstructed values for $l = 1, 2, 3$ , and expressed as

{\hat{W}}_{i + \frac{1}{2}}^{1} = p_{1} (x_{i + \frac{1}{2}}),

(30)

{\hat{W}}_{i + \frac{1}{2}}^{2} = \frac{1}{ξ_{2, 2}} p_{2} (x_{i + \frac{1}{2}}) - \frac{ξ_{1, 2}}{ξ_{2, 2}} {\hat{W}}_{i + \frac{1}{2}}^{1},

(31)

{\hat{W}}_{i + \frac{1}{2}}^{3} = \frac{1}{ξ_{3, 3}} p_{3} (x_{i + \frac{1}{2}}) - \frac{ξ_{1, 3}}{ξ_{3, 3}} {\hat{W}}_{i + \frac{1}{2}}^{1} - \frac{ξ_{2, 3}}{ξ_{3, 3}} {\hat{W}}_{i + \frac{1}{2}}^{2},

(32)

with

p_{1} (x_{i + \frac{1}{2}}) = {\bar{W}}_{i},

(33)

p_{2} (x_{i + \frac{1}{2}}) = \frac{- 1}{6} {\bar{W}}_{i - 1} + \frac{5}{6} {\bar{W}}_{i} + \frac{1}{3} {\bar{W}}_{i + 1},

(34)

\begin{matrix} p_{3} (x_{i + \frac{1}{2}}) = \frac{1}{30} {\bar{W}}_{i - 2} - \frac{13}{60} {\bar{W}}_{i - 1} + \frac{47}{60} {\bar{W}}_{i} \\ + \frac{9}{20} {\bar{W}}_{i + 1} + \frac{1}{20} {\bar{W}}_{i + 2} . \end{matrix}

(35)

Where $ξ_{1, 2} + ξ_{2, 2} = 1$ , $ξ_{1, 3} + ξ_{2, 3} + ξ_{3, 3} = 1$ , and $ξ_{2, 2}$ , $ξ_{3, 3} \neq 0$ . The linear weights are denoted by $ξ_{1, 2}, ξ_{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, $ξ_{1, 2} = 1 / 11$ , $ξ_{2, 2} = 10 / 11$ , $ξ_{1, 3} = 1 / 111$ , $ξ_{2, 3} = 10 / 111$ , and $ξ_{3, 3} = 100 / 111$ for the fifth order approximation, that is explained in Zhu and Shu.¹⁷

While, non linear weights are $ω_{l}$ given in equation (28) are expressed as follows

ω_{l} = \frac{{\bar{ω}}_{l}}{\sum_{l = 1}^{m} {\bar{ω}}_{l}}, {\bar{ω}}_{l} = ξ_{l, m} (1 + \frac{τ}{ε + B_{l}}), l = 1, . . ., m; m = 3 .

(36)

Here, $ξ_{l}$ represent the linear weights and $B_{l}$ represent smoothness indicators and in all the simulations $ε$ is taken as $10^{- 10}$ . While expression for given smoothness indicators in general form is expressed as follows

B_{m} = \sum_{j = 1}^{k} \int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} h^{2 j - 1} {(\frac{d^{j} p_{m} (x)}{d x^{j}})}^{2} dx, m = 2, 3,

(37)

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.¹⁷ To be more accurate, these smooth indicators are expressed as follow

B_{2} = \frac{13}{12} {({\bar{W}}_{i - 1} - 2 {\bar{W}}_{i} + {\bar{W}}_{i + 1})}^{2} + \frac{1}{4} {({\bar{W}}_{i - 1} - {\bar{W}}_{i + 1})}^{2},

(38)

\begin{matrix} B_{3} = {(η_{1} + \frac{1}{10} η_{3})}^{2} + \frac{13}{3} {(η_{2} + \frac{123}{455} η_{4})}^{2} \\ + \frac{781}{20} {(η_{3})}^{2} + \frac{1421461}{2275} (η_{4}), \end{matrix}

(39)

where

η_{1} = \frac{1}{120} (11 {\bar{W}}_{i - 2} - 82 {\bar{W}}_{i - 1} + 82 {\bar{W}}_{i + 1} - 11 {\bar{W}}_{i + 2}),

(40)

η_{2} = \frac{1}{56} (- 3 {\bar{W}}_{i - 2} + 40 {\bar{W}}_{i - 1} - 74 {\bar{W}}_{i + 1} + 40 {\bar{W}}_{i + 1} - 3 {\bar{W}}_{i + 2}),

(41)

η_{3} = \frac{1}{12} (- {\bar{W}}_{i - 2} + 2 {\bar{W}}_{i - 1} - 2 {\bar{W}}_{i + 1} + {\bar{W}}_{i + 2}),

(42)

η_{4} = \frac{1}{24} ({\bar{W}}_{i - 2} - 4 {\bar{W}}_{i - 1} + 6 {\bar{W}}_{i} - 4 {\bar{W}}_{i + 1} + {\bar{W}}_{i + 2}) .

(43)

To resolve the discontinuities efficiently, $B_{1}$ is expressed as follows

δ_{0} = {({\bar{W}}_{i} - {\bar{W}}_{i - 1})}^{2}, δ_{1} = {({\bar{W}}_{i + 1} - {\bar{W}}_{i},)}^{2}

(44)

{\bar{ξ}}_{0, 1} = {\begin{matrix} 1, & δ_{0} > δ_{1}, \\ 10, & otherwise \end{matrix}, {\bar{ξ}}_{1, 1} = 1 - {\bar{ξ}}_{0, 1},

(45)

ξ_{0, 1} = \frac{{\bar{ξ}}_{0, 1}}{{\bar{ξ}}_{0, 1} + {\bar{ξ}}_{1, 1}}, ξ_{1, 1} = 1 - ξ_{0, 1},

(46)

\begin{matrix} σ_{0} = ξ_{0, 1} (1 + \frac{{| δ_{0} - δ_{1} |}^{2}}{δ_{0} + ε}), σ_{1} = ξ_{1, 1} (1 + \frac{{| δ_{0} - δ_{1} |}^{2}}{δ_{1} + ε}), \\ σ = σ_{0} + σ_{1} . \end{matrix}

(47)

In equation (47), To avoid the denominator to become zero $ε$ is taken as a small positive number. Finally, we obtained

B_{1} = \frac{1}{σ^{2}} {(σ_{0} ({\bar{W}}_{i} - {\bar{W}}_{i - 1}) + σ_{1} ({\bar{W}}_{i + 1} - {\bar{W}}_{i}))}^{2} .

(48)

The term $τ$ in equation (36) is stated as

τ = (\frac{\sum_{l = 1}^{2} | B_{3} - B_{l} |}{2}) .

(49)

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 $\int_{x_{i - \frac{1}{2}}}^{x_{i + \frac{1}{2}}} S (W) dx$ is denoted by $S_{i}$ and given as follows

S_{i} = {(0, 0, \frac{32}{D^{2}} ({(u_{mix})}_{j + \frac{1}{2}}) (μ_{mix})_{j + \frac{1}{2}}) + {(u_{mix})}_{j - \frac{1}{2}}) (μ_{mix})_{j - \frac{1}{2}})))}^{t},

(50)

where $t$ denotes the transpose. Thus, we left the system of ODEs which are

\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},

(51)

the above system can also be rewritten as

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

(52)

Now, the third order TVD RK method¹⁷ is used to compute the system equation (52), as follow

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

(53)

where, by using the relation $\frac{CFL * dx}{max (| λ_{1} |, | λ_{2} |, | λ_{3} |)}$ 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³³ 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.³⁴ 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

(ρ_{g}, ρ_{l}, u, a_{g}, a_{l}) (x, 0) = {\begin{matrix} (1, 1, 0, 500, 100), & x > 0, \\ (6, 1, 0, 500, 100), & x < 0 . \end{matrix}

(54)

For test problem 1, $p$ , $α$ the pressure and speed of sound are determined by using the following relations

\sqrt{\frac{p}{ρ_{g} a_{g} + ρ_{l} a_{l}}} and a_{g}^{2} ρ_{g} + a_{l}^{2} ρ_{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.¹⁰ Moreover, 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) = {\begin{matrix} (0.55, 10.37, 80450), & x < 50, \\ (0.55, 0.561, 24282), & x > 50 . \end{matrix}

Figure 1.

The solution profiles of test problem 1 at $t = 0.5 s$ .

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.

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) = {\begin{matrix} (0.35, 1.868, 192170), & x < 50, \\ (0.30, 14.47, 196690), & x > 50 . \end{matrix}

Figure 2.

The solution profiles of test problem 2 at $t = 1.0 s$ .

Table 2.

$L^{1}$ -errors for Test problem 2.

No. of cells	Gas velocity	Gas velocity	Gas velocity
	$L^{1}$ -error (simple WENO scheme)	$L^{1}$ -error (MCUP scheme)	$L^{1}$ -error (MR-WENO-AUSMV scheme)
$100$	$1.47$	$2.21$	$1.21$
$200$	$0.91$	$1.41$	0.81
$400$	$0.15$	$0.93$	0.10
$800$	$0.087$	$0.13$	0.056
$1600$	$0.011$	$0.079$	0.009

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 are resolved by derived numerical scheme (MR WENO) more accurately and efficiently then the high-resolution numerical scheme (MCUP).

Test problem 4: we have taken this problem from the article.⁶ This 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.4 s$ and $t = 1.3 s$ are shown in top row of the Figure 4. While, the profiled solution for the pressure at times $t = 0.4 s$ and $t = 1.3 s$ 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.⁶ The source term in equation (2) and the gas-slip relationship in equation (7) in the the algebraic form are used by taking values of $u_{d} = 0.0 m / s$ and $C_{o} = 1.0$ . Initially, we have considere $10^{5}$ Pa pressure in the horizontal pipe whose length is $1000 m$ having diameter $0.1 m$ . 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 = 250 m$ , $x = 750 m$ , respectively, and the bottom row of figures show gas mass flow rates profiles at different position $x = 250 m$ , $x = 750 m$ within the solution space. Our final results are in close agreement with those presented in Fjelde and Karlsen,⁶ 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.⁶ The 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.³⁵ In 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.³⁵ th 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

u_{d} = 0.35 \sqrt{\frac{g (ρ_{l} - ρ_{g}) D}{ρ_{l}}} .

Figure 3.

The solution profiles of test problem 3 at $t = 1.0 s$ .

Figure 4.

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

Figure 5.

The solution profiles of test problem 5.

This type of drift velocity³⁵ 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. Therefore the gas volume fraction and gas phase velocities increases as shown in the Figure 6.

Figure 6.

Solution profiles obtained at time $t = 250 s$ .

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 fifth-order 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.

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 iD

Asad Rehman

References

Pelanti

LeVeque

RJ.

High-resolution finite volume methods for dusty gas jets and plumes. SIAM J Sci Comput 2006; 28: 1335–1360.

Dobran

Neri

Macedonio

Numerical simulation of collapsing volcanic columns. J Geophys Res 1993; 98: 4231–4259.

Dufek

Bergantz

GU.

Controls of basal particle loss on pyroclastic flow transport. Geochem Geophy Geosys 2007; 8: 12. DOI: 10.1029/2007GC001722.

Zuber

Findlay

JA.

Average volumetric concentration in two-phase flow systems. J Heat Transf 1965; 87: 453–468.

Romate

JE.

An approximate Riemann solver for a two-phase flow model with numerically given slip relation. Comput Fluids 1998; 27: 455–477.

Fjelde

Karlsen

KH.

High-resolution hybrid primitive-conservative upwind schemes for the drift-flux model. Comput Fluids 2002; 31: 335–367.

Steinar

Fjelde

KK.

Hybrid flux-splitting schemes for a two-phase flow model. J Comput Phys 2002; 175: 674–701.

Baudin

Berthon

Coquel

, et al. A relaxation method for two-phase flow models with hydrodynamic closure law. Numer Math 2005; 99: 411–440.

Munkejord

Evje

Flatten

The multi-stage centered-scheme approach applied to a drift-flux two-phase flow model. Int J Numer Methods Fluids 2006; 52: 679–705.

10.

Evje

Fjelde

KK.

Relaxation schemes for the calculation of two-phase flow in pipes. Math Comput Model 2002; 36: 535–567.

11.

Evje

Fjelde

KK.

On a rough AUSM scheme for a one-dimensional two-phase model. Comput Fluids 2003; 32: 1497–1530.

12.

Steinar

Tore

Munkejord

ST.

A WIMF scheme for the drift-flux two-phase flow model. Citeseer, 2006.

13.

Danish

Ashiq

Asad

, et al. The space-time conservation element and solution element scheme for simulating two-phase flow in pipes. Adv Mech Eng 2019; 11: 168781401989835.

14.

Paula

Valdes

JP.

Analysis of the drift flux in two-phase gas-liquid slug-flow along horizontal and inclined pipelines. Int J Heat Mass Transf 2018; 5: 53–70.

15.

Rehman

Meraj

MA.

A limiter free adaptive THINC-BVD scheme with AUSMV numerical flux for solving the drift flux model. Int J Non Linear Mech 2022; 141: 103931.

16.

Welahettige

Chiristian

Bernt

Development of central-upwind-WENO scheme for two-phase 1-D drift-flux model. Scand Simul Soc 2022: 345–352.

17.

Zhu

Shu

C-W.

A new type of multi-resolution WENO schemes with increasingly higher order of accuracy. J Comput Phys 2018; 315: 659–683.

18.

Evje

Fjelde

KK.

Hybrid flux-splitting schemes for a two-phase flow model. J Comput Phys 2002; 175: 674–701.

19.

Liu

Osher

Chan

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

20.

Harten

Osher

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

21.

Shu

C-W

Osher

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

22.

Shu

C-W

Osher

Efficient implementation of essentially non-oscillatory shock capturing schemes, II. J Comput Phys 1989; 83: 32–78.

23.

Shu

C-W.

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

24.

Balsara

Shu

C-W.

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

25.

Shu

C-W.

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

26.

Xing

Shu

C-W.

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

27.

Shu

C-W.

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

28.

Zhang

Shu

C-W.

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

29.

Zhu

Qiu

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

30.

Zhu

Shu

C-W.

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

31.

Liou

Steffen

Jr.

A new flux splitting scheme. J Comput Phys 1993; 107: 23–39.

32.

Wada

Liou

MS.

An accurate and robust flux splitting scheme for shock and contact discontinuities. SIAM J Sci Comput 1997; 18: 633–657.

33.

Kurganov

Lin

CT.

On the reduction of numerical dissipation in central-upwind schemes, Commun Comput Phys 2007; 2: 141–163.

34.

Kuila

Raja Sekhar

Zeidan

A robust and accurate Riemann solver for a compressible two-phase flow model. App Math Comput 2015; 265: 681–695.

35.

Velmurugan

Dynamic simulation of gas expansion in Marine Risers. Master’s Thesis, NTNU, Norway, 2016.