The space-time conservation element and solution element scheme for simulating two-phase flow in pipes

Abstract

In this article, the space-time conservation element and solution element scheme is extended to simulate the unsteady compressible two-phase flow in pipes. The model is non-conservative and the governing equations consist of three equations, namely, two mass conservation equations for each phase and one mixture-momentum equation. In the third equation, the non-conservative source term appears, which describes the sum of gravitational and frictional forces. The presence of source term and two mass conservation equations in considered model offers difficulties in developing the accurate and robust numerical techniques. The suggested space-time conservation element and solution element numerical scheme resolves the volume-contact discontinuities efficiently. Furthermore, the modified central upwind scheme is also extended to solve the same two-phase flow model. The number of test problems is considered, and the results obtained by space-time conservation element and solution element scheme are compared with the solutions of modified central upwind scheme. The numerical results show better performance of the space-time conservation element and solution element method as compare to the modified central upwind scheme.

Keywords

Drift-flux model space-time conservation element and solution element scheme source term shock wave volume-fraction contact discontinuity

Introduction

Two-phase flows are widely observed in natural environment, such as snowy or rainy winds, typhoons, water and air pollution, and volcanic eruptions. They are also extensively encountered in nuclear power plants, combustion engines, bio-medical engineering, food processing industry, and many more.^1–3 Because of these wide range applications, many researchers had developed several two-phase models and numerical methods to study the dynamics of such flows. These two-phase flow models inherit several numerical difficulties since each phase is considered separately and the model comprises two sets of conservation of mass, momentum, and energy. For example, the presence of two momentum equations causes difficulties such as loss of hyperbolicity and uncertainties in specifying interfacial interaction terms between the two phases. Here, we are interested in the model that elaborates the compressible two-phase flow in gas and liquid horizontal pipelines. However, above-mentioned difficulties in two-phase flow models are remarkably reduced by formulating the drift-flux model, in which the mixture momentum equation is utilized to describe the motion of whole mixture. This model is a powerful tool for gaining insight into the flow processes where oil and gas are transported simultaneously out of a reservoir. Without modeling of such flows in oil industry, it is very difficult to describe effects that may arise within a well bore over time to weaken or improve production.

In this article, we consider the drift-flux model⁴ that consists of two mass conservation equations for each phase and one momentum equation for the mixture. The drift-flux model was initially designed by Zuber and Findlay⁵ and has been improved by many researchers.^4,6–⁹ The considered drift-flux model is derived from the two-fluid model by adding the two momentum equations for each phase.¹⁰ Now, more difficult terms associated with phase interaction are canceled out, and the missing information is obtained using the kinematic constitutive equation,⁷ which relates the phase velocities at any point. Still, the source term that is sum of gravitational and frictional forces is presented in the considered model. This model explains the isothermal liquid-gas flow in the long pipeline where the flow behavior is averaged and perpendicular to the pipe axis. Hence, the resulting model is one-dimensional (1D) in the direction of axis.

Due to importance of the drift-flux model, many numerical schemes had been designed and extended to investigate this model in the literature. First of all, Romate⁶ has developed an approximate Roe-type Riemann solver for computing the drift-flux model. Afterward, hybrid flux-splitting and relaxation-type numerical schemes are designed to solve the same model in Evje and Fjelde,^7,11 respectively. Meanwhile, Fjelde and Karlsen⁴ have developed high-resolution hybrid upwind scheme to analyze the considered model. In this scheme, for a smooth region, a simple non-conservative Monotonic Upwind Scheme for Conservation Laws (MUSCL) scheme and, in a region of strong discontinuities, high-resolution conservative scheme are utilized. Next, in Evje and Fjelde,¹² another splitting method, namely, advection upstream splitting method, is extended to analyze the drift-flux model.⁷ Furthermore, for solving the same model, semi-implicit relaxation scheme is proposed in Baudin et al.⁸ and a numerical scheme weakly implicit mixture flux (WIMF) is extended in Steinar et al.¹³ Subsequently, the multi-stage approach (MUSTA) is used to develop the centered numerical scheme for investigating the drift-flux model in Munkejord et al.⁹ Recently, Paula and Valdes¹⁴ analyze the drift-flux model in two-phase slug-flow in horizontal and inclined pipelines using experimental non-Newtonian and Newtonian approaches.

In this article, the space-time conservation element and solution element (CE/SE) scheme¹⁵ is extended for simulating the compressible two-phase flow in horizontal oil and gas pipelines. This scheme is entirely different from the schemes which have been applied before to solve the considered model. This method has many distinct features, such as treatment of space and time at the same step, introduction of conservation elements (CEs) and solution elements (SEs), shock capturing approach without utilizing Riemann solvers, and the use of staggered grid. Besides these distinct features, the suggested scheme is distinguished by the simplicity of its conceptual basis-conservation of flux in time and space; for details, see Chang.¹⁵ Various applications of CE/SE scheme in different areas affirm the scheme’s generality, robustness, and effectiveness.^15–24 Later on, H Shen and colleagues^25,26 designed a new upwind CE/SE scheme which is based on the “a” scheme (the original CE/SE scheme). In this new scheme, the numerical dissipation is added through the upwind procedure. This upwind CE/SE scheme preserves almost all features of the original CE/SE scheme. For extensive detail, the reader is referred to the literature.^27–31 The number of test problems is considered to show that the suggested scheme is highly robust, gives better resolutions of the sharp volume-fraction contact discontinuities, and preserves the positivity of flow variables such as fluid densities, volume fractions, and pressure. For checking the accuracy of proposed numerical scheme, we have extended the modified central upwind scheme (CUP)³² for solving the considered drift-flux model, and the results obtained from CE/SE scheme are compared with those of modified CUP.

The rest of article is organized as follows. In section “Drift-flux model,” the drift-flux model is given and data-dependent terms are described. The CE/SE scheme for drift-flux model is described in section “Construction of CE/SE scheme for 1D drift-flux model.” In section “Numerical test problems,” different interesting test problems are included, to validate and compare the results of suggested numerical schemes. Finally, a conclusion of this article is drawn in section “Conclusion.”

Drift-flux model

In this section, we present the mathematical form and eigen-structure of the drift-flux model. Also, we describe the submodels that are used in this article. We will restrict ourselves to consider liquid-gas isothermal flow in the horizontal pipeline. The fundamental equations for the considered drift-flux model are written as

\begin{matrix} \partial_{t} (a_{l} ρ_{l}) + \partial_{x} (a_{l} ρ_{l} u_{l}) = 0 \\ \partial_{t} (a_{g} ρ_{g}) + \partial_{x} (a_{g} ρ_{g} u_{g}) = 0 \\ \partial_{t} (a_{l} ρ_{l} u_{l} + a_{g} ρ_{g} u_{g}) + \partial_{x} (a_{l} ρ_{l} u_{l}^{2} + a_{g} ρ_{g} u_{g}^{2} + p) = - q \end{matrix}

(1)

Here, $a_{l}$ and $a_{g}$ are volume fractions, $ρ_{l}$ and $ρ_{g}$ are densities, and $u_{l}$ and $u_{g}$ are velocities. Subscripts $g$ and $l$ denote the gas and liquid, respectively. The mutual pressure for gas and liquid is denoted by $p$ . In third equation, $q$ on the right-hand side is the source term and defined as

q = F_{w} + F_{g} = \frac{32 u_{mix} μ_{mix}}{d^{2}} + gsin φ (a_{l} ρ_{l} + a_{g} ρ_{g})

(2)

Here, $d$ denotes the inner diameter of pipeline, $φ$ represents the inclination, and g is the gravitational constant. In this study, we have considered horizontal pipes only, thus $φ = 0$ ; however, there are some experimental and theoretical studies available in the literature which study liquid-gas flow in the inclined pipelines.³³ The mixture quantities in frictional forces $F_{w}$ , average velocity $u_{mix}$ , and viscosity $μ_{mix}$ are defined as

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

(3)

The considered model contains seven unknowns $u_{l}, u_{g}$ , $ρ_{l}, ρ_{g}$ , $a_{l}, a_{g}$ , $p$ and three equations. Thus, to obtain the solutions of the system, additional constraints are needed. For this purpose, we have the following relations:

1. Volume fractions are related by the following relation

a_{l} + a_{g} = 1

(4)

2. We use submodel for density of liquid as

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

(5)

and submodel for density of gas as follows

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

(6)

Here, $α_{l} and α_{g}$ are, respectively, sound velocities in fluid phase and $ρ_{l, 0}$ is reference liquid density with corresponding reference liquid pressure $p_{l, 0}$ .

3. We use the following algebraic form of gas slip relation for computational purposes

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

(7)

Here, $C_{o}$ is the flow dependent parameter, and $u_{d}$ and $u_{m}$ denote drift-gas and drift-mixture velocity, respectively. The values of constant parameters used in equations (3)–(7) are given in Table 1.

Table 1.

Parameters used in numerical computation of the drift-flux model.

Parameters	Value	Description
$μ_{l}$	$5 \times 10^{- 2} Pa s$	Liquid viscosity
$μ_{g}$	$5 \times 10^{- 6} Pa s$	Gas viscosity
$α_{l}$	$1000 m / s$	Sound velocity in liquid phase
$α_{g}$	$316.22 m / s$	Sound velocity in gas phase
$ρ_{l, 0}$	$1000 kg / m^{3}$	Reference density liquid
$p_{l, 0}$	$10^{5} Pa$	Reference pressure liquid

Eigen-structure

The system of equation (1) can be written as

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

(8)

where

\begin{matrix} 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}), Q (W) = (\begin{matrix} 0 \\ 0 \\ - q \end{matrix}) \end{matrix}

with $W = (w_{1}, w_{2}, w_{3}) = (a_{l} ρ_{l}, a_{g} ρ_{g}, a_{l} ρ_{l} u_{l} + a_{g} ρ_{g} u_{g})$ , $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 $Q (W) = (q_{1}, q_{2}, q_{3}) = (0, 0, - q)$ . Now, equation (8) can be rewritten as

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

(9)

Note that the passive variable pressure can be obtained from the conservative variables $w_{1}$ and $w_{2}$ . Clearly, equation (9) shows that the components of physical flux vector cannot be written completely in terms of conserved variables. Hence, we use some assumptions as described in Fjelde and Karlsen⁴ and obtain the approximate sound velocity. For more comprehensive detail about the properties of drift-flux model, reader is referred to the literature.^4,7 Here, we present only main assumptions that are directly involved to obtain the approximate sound velocity.

First, by assuming $C_{o} = 1$ and $u_{d} = 0$ in equation (7), we get the no-slip condition, that is, $u_{g} = u_{l} = u$ . Now, using this assumption, the common velocity $u$ can be defined in terms of the conserved variable as $u (W) = w_{3} / w_{1} + w_{2}$ . Using this expression of $u (W)$ , equation (9) becomes

\partial_{t} (\begin{matrix} w_{1} \\ w_{2} \\ w_{3} \end{matrix}) + \partial_{x} (\begin{matrix} \frac{w_{1} w_{3}}{w_{1} + w_{2}} \\ \frac{w_{2} w_{3}}{w_{1} + w_{2}} \\ \frac{w_{1} w_{3}^{2}}{{(w_{1} w_{2})}^{2}} + \frac{w_{2} w_{3}^{2}}{{(w_{1} + w_{2})}^{2}} + p (w_{1}, w_{2}) \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ - q \end{matrix})

(10)

The above system is the required expression for finding the Jacobian matrix. Now, by applying the definition of $w_{1}$ and $w_{2}$ , equation (4) in terms of $ρ$ and $w$ can be written as

\frac{w_{1}}{ρ_{l} (p)} + \frac{w_{2}}{ρ_{g} (p)} = 1

(11)

Using the values of phase densities $ρ_{l, g} (p)$ as defined in equations (5) and (6), equation (11) becomes

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

(12)

By rearranging equation (12) with respect to pressure $p$ , we get a second-order polynomial of the type

p^{2} + Bp + C = 0

(13)

where $B = α_{l}^{2} (ρ_{l, 0} - (p_{l, o} / α_{l}^{2}) - w_{1} - (α_{g} / α_{l})^{2} w_{2})$ and $C = w_{2} (α_{l} α_{g})^{2} ((p_{l, o} / α_{l}^{2}) - ρ_{l, 0})$ . The negative root yields negative pressures for all $w_{1}$ and $w_{2}$ for solving pressure, so here only positive root will be considered for the physical solutions. The pressure is then given by

p (w_{1}, w_{2}) = \frac{- B + \sqrt{B^{2} - 4 C}}{2}

(14)

Now, from the flux vector in equation (10), we obtain the Jacobian matrix as follows

J = (\begin{matrix} \frac{w_{2}}{w_{1} + w_{2}} u - \frac{w_{1}}{w_{1} + w_{2}} u \frac{w_{1}}{w_{1} + w_{2}} \\ - \frac{w_{2}}{w_{1} + w_{2}} u \frac{w_{1}}{w_{+} w_{2}} u \frac{w_{2}}{w_{1} + w_{2}} \\ - u^{2} + \frac{\partial p}{\partial w_{1}} - u^{2} + \frac{\partial p}{\partial w_{2}} 2 u \end{matrix})

(15)

with $u = w_{3} / (w_{1} + w_{2})$ .

Second, we assume that $a_{g} ρ_{g} << a_{l} ρ_{l}$ , that is, $w_{2} << w_{1}$ . Now, the Jacobian matrix (15) becomes

J = (\begin{matrix} 0 & - u & 1 \\ 0 & u & 0 \\ - u^{2} + \frac{\partial p}{\partial w_{1}} & - u^{2} + \frac{\partial p}{\partial w_{2}} & 2 u \end{matrix})

The corresponding eigenvalues are

λ_{1} = u - α, λ_{2} = u, λ_{3} = u + α

(16)

where $λ_{2}$ represents the velocity of the moving fluid while $λ_{1}$ and $λ_{3}$ represent the velocity of pressure pulses traveling in opposite directions. Thus, approximate sound velocity $α$ is obtained with the assumptions of $a_{g} ρ_{g} << a_{l} ρ_{l}$ and $ρ_{l} \approx constant$ ; for details, see Fjelde and Karlsen⁴

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

(17)

Construction of CE/SE scheme for 1D drift-flux model

The CE/SE numerical has entirely different concept and methodology from the well-established numerical techniques such as the finite difference, finite volume, and finite element methods. In this section, the 1D CE/SE scheme¹⁵ is constructed for 1D drift-flux model (1). For the construction of CE/SE scheme, we rewrite equation (8) in component form as

\frac{\partial w_{m}}{\partial t} + \frac{\partial f_{m}}{\partial x} - Q_{m} = 0, m = 1, 2, 3

(18)

Let $x_{0} = t$ , $x_{1} = x$ be the coordinates of a two-dimensional Euclidean space $(E_{2})$ . By applying the Gauss divergence theorem, equation (18) is equivalent to the integral equation

\oint_{S (V)} W \cdot d S = \int_{V} Q_{m} dV

(19)

where $W = (w_{m}, f_{m})^{T}$ , $m = 1, 2, 3,$ that is, for each $m = 1, 2, 3$ , the components of the vector $W$ in the $t$ and $x$ directions are $w_{m}$ and $f_{m}$ , respectively, and $d S \overset{def}{=} d σ n$ with $d σ$ and $n$ being the area and unit outward normal vector of a surface element on $S (V)$ . Now, equation (19) is enforced over a space-time domain, called CE, that allows the discontinuities of flow variables. The actual numerical integration is carried out in a discrete manner using SEs. An SE is a different space-time region in which the flow variables are assumed to be smooth, thus discretization of the flow variables with a prescribed order of accuracy can be performed.

We denote the computational domain by $Ω$ as the set of mesh points $(j, n)$ in $E_{2}$ space where $n = 0, 1 / 2, 1, 3 / 2, \dots$ and for each $n$ , $j = 0, \pm 1 / 2, \pm 1, \pm 3 / 2, \dots$ . For each mesh point $(j, n)$ , there is an SE associated with it—shown in Figure 1—as the interior of the space-time region denoted by a dashed curve. It consists of a horizontal line segment, a vertical line segment, and their immediate neighborhood. The exact size of the neighborhood does not matter. For any $(x, t) \in SE (j, n), w_{m} (x, t), f_{m} (x, t)$ , and $W (x, t)$ are approximated by $w_{m}^{*} (x, t; j, n), f_{m}^{*} (x, t; j, n)$ , and $W^{*} (x, t; j, n)$ , respectively, as follows

w_{m}^{*} (x, t; j, n) = {(w_{m})}_{j}^{n} + {(w_{mt})}_{j}^{n} (t - t^{n}) + {(w_{mx})}_{j}^{n} (x - x_{j})

(20)

and

f_{m}^{*} (x, t; j, n) = {(f_{m})}_{j}^{n} + {(f_{mt})}_{j}^{n} (t - t^{n}) + {(f_{mx})}_{j}^{n} (x - x_{j})

(21)

Figure 1.

Space-time staggered grid near SE (j, n) (a) Space-time staggered grid near SE(j,n), (b) CE_(-)(j,n) and CE_(+)(j,n), and (c) CE(j,n).

Moreover

Q_{m}^{*} (x, t; j, n) = Q_{m} (w_{m}^{*}, w_{mx}^{*})

(22)

By employing chain rule, we obtain

{(f_{mx})}_{j}^{n} \overset{def}{=} \sum_{k = 1}^{3} {(f_{m, k})}_{j}^{n} {(w_{kx})}_{j}^{n}

(23)

{(f_{mt})}_{j}^{n} = \sum_{k = 1}^{3} {(f_{m, k})}_{j}^{n} {(w_{kt})}_{j}^{n}

(24)

where

f_{m, k} = \frac{\partial f_{m}}{\partial w_{k}}, m, k = 1, 2, 3

(25)

The Jacobian matrix is formed by $f_{m, k}, m, k = 1, 2, 3$ , with $m$ and $k$ being the row and column indices, respectively. Note that $(w_{m})_{j}^{n}$ , $(w_{mx})_{j}^{n}$ , and $(w_{mt})_{j}^{n}$ are constant in $SE (j, n)$ . These are numerical analogues of the values of $w_{m}$ , $\partial w_{m} / \partial x$ , and $\partial w_{m} / \partial t$ at $(x_{j}, t^{n})$ , respectively. Since $W = (w_{m}, f_{m})^{T}$ , one can write

W^{*} (x, t : j, n) = {(w_{m}^{*} (x, t : j, n), f_{m}^{*} (x, t : j, n))}^{T}

(26)

Moreover, due to the smoothness assumption of variables for any $(x, t) \in SE (j, n), w_{m}^{*} (x, t; j, n)$ , $f_{m}^{*} (x, t; j, n)$ , and $Q_{m}^{*} (x, t : j, n)$ satisfy equation (1), that is

\frac{\partial w_{m}^{*} (x, t; j, n)}{\partial t} + \frac{\partial f_{m}^{*} (x, t; j, n)}{\partial x} = Q_{m}^{*} (x, t; j, n)

(27)

Using equations (20) and (21), equation (27) is equivalent to

{(w_{mt})}_{j}^{n} = - {(f_{mx})}_{j}^{n} + Q_{m} ({(w_{m})}_{j}^{n}, {(w_{mx})}_{j}^{n})

(28)

We notice that $(f_{m})_{j}^{n}$ is function of $(w_{m})_{j}^{n}$ , $(f_{mx})_{j}^{n}$ is function of $(w_{m})_{j}^{n}$ and $(w_{mx})_{j}^{n}$ , and $(f_{mt})_{j}^{n}$ is function of $(w_{m})_{j}^{n}$ and $(w_{mt})_{j}^{n}$ . As a result, $(w_{m})_{j}^{n}$ and $(w_{mx})_{j}^{n}$ are the only independent discrete variables to be calculated in the current marching scheme.

Let $E_{2}$ be divided into non-overlapping rectangular regions (see Figure 1(a)) referred to as CEs. As depicted in Figure 1(b), two CEs, that is, $C E_{-} (j, n)$ and $C E_{+} (j, n)$ , are associated with each interior mesh point $(j, n) \in Ω$ . These CEs are referred to as basic conservation elements (BCEs). Contrarily, $CE (j, n)$ in Figure 1(c), which is the union of $C E_{-} (j, n)$ and $C E_{+} (j, n)$ , is called compounded conservation element (CCE).

Note that, among the line segments forming the boundary of $C E_{-} (j, n)$ , $\bar{AB}$ and $\bar{AD}$ belong to $SE (j, n)$ , while $\bar{CB}$ and $\bar{CD}$ belong to $SE (j - 1 / 2, n - 1 / 2)$ . Similarly, the boundary of $C E_{+} (j, n)$ belongs to either $SE (j, n)$ or $SE (j + 1 / 2, n - 1 / 2)$ . As a result, by imposing two conservation conditions at each $(j, n) \in Ω$ , we obtain

\oint_{S (C E_{\pm} (j, n))} W^{*} \cdot d S = \oint_{C E_{\pm} (j, n)} Q_{m} dV, m = 1, 2, 3

(29)

According to this equation, the total flux leaving boundary of any BCE is zero. Because the surface integration over any interface separating two neighboring BCEs is evaluated using the information from a single SE, obviously, the local conservation relation (29) leads to a global flux conservation relation, that is, the total flux of $W^{*}$ leaving the boundary of any space-time region that is the union of any conservation of BCEs will also vanish. In particular, because $CE (j, n)$ is the union of $C E_{-} (j, n)$ and $C E_{+} (j, n)$

\int_{S (CE (j, n))} W^{*} \cdot d S = \int_{CE (j, n)} Q_{m} dV, m = 1, 2, 3, (j, n) \in Ω

(30)

must follow from equation (29).

Using equation (29) along with equations (20)–(22) and (26), one obtains

\begin{matrix} {(w_{m})}_{j}^{n} = \frac{1}{2} [{(w_{m})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} + {(w_{m})}_{j - \frac{1}{2}}^{n - \frac{1}{2}}] \\ - \frac{Δ x}{8} [{(w_{mx})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} - {(w_{mx})}_{j - \frac{1}{2}}^{n - \frac{1}{2}}] - \frac{Δ t}{2 Δ x} [{(f_{m})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} - {(f_{m})}_{j - \frac{1}{2}}^{n - \frac{1}{2}}] \\ - \frac{Δ t^{2}}{8 Δ x} [{(f_{mt})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} - {(f_{mt})}_{j - \frac{1}{2}}^{n - \frac{1}{2}}] + \frac{Δ t}{4} [{(Q_{m})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} + {(Q_{m})}_{j - \frac{1}{2}}^{n - \frac{1}{2}}] \end{matrix}

(31)

or above equation can be written as

{(w_{m})}_{j}^{n} = \frac{1}{2} [{(w_{m})}_{j - \frac{1}{2}}^{n - \frac{1}{2}} + {(w_{m})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} + {(S_{m})}_{j - \frac{1}{2}}^{n - \frac{1}{2}} - {(S_{m})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} + \frac{Δ t}{2} ({(Q_{m})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} + {(Q_{m})}_{j - \frac{1}{2}}^{n - \frac{1}{2}})]

(32)

where

{(S_{m})}_{j \pm \frac{1}{2}}^{n - \frac{1}{2}} = \frac{Δ x}{4} {(w_{mx})}_{j \pm \frac{1}{2}}^{n - \frac{1}{2}} + \frac{Δ t}{Δ x} {(f_{m})}_{j \pm \frac{1}{2}}^{n - \frac{1}{2}} + \frac{Δ t^{2}}{4 Δ x} {(f_{mt})}_{j \pm \frac{1}{2}}^{n - \frac{1}{2}}

(33)

Note that $q_{1} = 0 = q_{2} = 0$ and $q_{3} \neq 0$ . Moreover

{(q_{3})}_{j + \frac{1}{2}}^{n - \frac{1}{2}} + {(q_{3})}_{j - \frac{1}{2}}^{n - \frac{1}{2}} = \frac{32}{d^{2}} ({(h)}_{j - \frac{1}{2}}^{n - \frac{1}{2}} {(θ)}_{j - \frac{1}{2}}^{n - \frac{1}{2}} + {(h)}_{j + \frac{1}{2}}^{n - \frac{1}{2}} {(θ)}_{j + \frac{1}{2}}^{n - \frac{1}{2}})

(34)

where $h$ is used in place of $u_{mix}$ and $θ$ is used in place of $μ_{mix}$ , which are given in equation (2).

The numerical oscillations near a discontinuity can be suppressed using the following limiting formulations for the slopes of conservative variables

{(w_{mx})}_{j}^{n} = U_{m} ({(w_{m x_{-}})}_{j}^{n}, {(w_{m x_{+}})}_{j}^{n}; ζ), m = 1, 2, 3

(35)

Here, $ζ \geq 0$ an adjustable constant (usually $ζ = 1$ or $α = 2$ ). In this article, we set $ζ = 1$ and

U_{m} (x_{-}, x_{+}; ζ) = \frac{| x_{+} |^{ζ} x_{-} + | x_{-} |^{ζ} x_{+}}{| x_{+} |^{ζ} + | x_{-} |^{ζ}}

(36)

Moreover

{(w_{m x_{+}})}_{j}^{n} = \frac{{({w'}_{m})}_{j + \frac{1}{2}}^{n} - {(w_{m})}_{j}^{n}}{Δ x / 2}, {(w_{m x_{-}})}_{j}^{n} = \frac{{(w_{m})}_{j}^{n} - {({w'}_{m})}_{j - \frac{1}{2}}^{n}}{Δ x / 2}

(37)

and

{({w'}_{m})}_{j \pm \frac{1}{2}}^{n} = {(w_{m})}_{j \pm \frac{1}{2}}^{n - \frac{1}{2}} + \frac{Δ t}{2} {(w_{mt})}_{j \pm \frac{1}{2}}^{n - \frac{1}{2}}, l = 1, 2, 3

(38)

Equations (32) and (35)–(38) constitute the CE/SE solver for 1D drift-flux model. This completes the construction of proposed solver.

Numerical test problems

In this section, several test problems are presented for the drift-flux flow model. The obtained results of CE/SE scheme are also compared with the results of modified CUP.³² Furthermore, the reference solutions are obtained using the CE/SE numerical scheme on uniform $5000$ grid cells.

The first two test problems are taken from Kuila et al.³⁴ for checking the accuracy of proposed numerical schemes. Our proposed numerical schemes show better performance in resolving strong discontinuities as compared to the designed numerical schemes in Kuila et al.³⁴ The pressure $p$ and speed of sound $α$ are calculated using the relations $a_{g}^{2} ρ_{g} + a_{l}^{2} ρ_{l}$ and $\sqrt{p / (ρ_{g} a_{g} + ρ_{l} a_{l})}$ , respectively. Furthermore, the time step $dt$ is computed using the relation $CFL * dx / max (| λ_{1} |, | λ_{2} |, | λ_{3} |)$ , where $CFL$ denotes the Courant–Friedrichs–Lewy condition. In the first two test problems, we consider $CFL = 0.6$ . Furthermore, gas and liquid phase velocities are taken the same, that is, no-slip condition is used, as described in Kuila et al.³⁴

Test problem 1

The left and right states of the Riemann problem are

\begin{matrix} (ρ_{g, L} = 50, ρ_{l, L} = 100, a_{g, L} = 500.0, a_{l, L} = 100.0 u_{g, L} = - 100.0, u_{l, L} = - 100.0 m / s) \\ (ρ_{g, R} = 50, ρ_{l, R} = 1000, a_{g, R} = 500.0, a_{l, R} = 100.0 u_{g, R} = 100.0, u_{l, R} = 100.0 m / s) \end{matrix}

(39)

Here, the subscripts $R$ and $L$ represent the right and left states, respectively. In this test problem, $a_{g}$ and $a_{l}$ represent the compressibility factors of gas and liquid, respectively. The solution profiles of the problem computed on $200$ grid cells at time $t = 1.85 s$ are given in Figure 2, which shows that the left wave is a rarefaction wave and a contact discontinuity, and right wave is also a rarefaction wave.

Figure 2.

The solution of Riemann problem at $t = 0.185 s$ by CE/SE and CUP schemes.

Test problem 2

The left and right states of the Riemann problem are

\begin{matrix} (ρ_{g, L} = 0.1, ρ_{l, L} = 0.1, a_{g, L} = 700.0, a_{l, L} = 300.0 u_{g, L} = 200.0, u_{l, L} = 200.0 m / s) \\ (ρ_{g, R} = 0.1, ρ_{l, R} = 0.1, a_{g, R} = 700.0, a_{l, R} = 300.0 u_{g, R} = - 200.0, u_{l, R} = - 200.0 m / s) \end{matrix}

(40)

The solution profiles of test problem are given in Figure 3 at time $t = 0.015$ and obtained on 200 grid cells. Figure 3 shows that the left wave is a shock wave, central wave is a contact discontinuity, and right wave is also a shock wave.

Figure 3.

Numerical results of Riemann problem at $t = 0.015 s$ by CE/SE and CUP schemes.

The next two test problems are taken from Fjelde and Karlsen.⁴ In these problems, we consider the liquid density $ρ_{l} = 1000 kg / m^{3}$ and use the gas slip-condition (7), with $C_{o} = 1.07$ and $u_{d} = 0.216 m / s$ . The density of gas is obtained using relation (6), and again for these problems, we consider $F_{g} = 0$ . Furthermore, the time step $dt$ is computed using the relation $CFL * dx / max (| λ_{1} |, | λ_{2} |, | λ_{3} |)$ , where the value of $α$ is taken from equation (17) and set $CFL = 0.4$ . In both test problems, the length of horizontal pipe is taken $100 m$ and solution profiles are computed on 400 grid cells. Initially, the pipe is separated in a left and right state at $x_{0} = 50 m$ .

Test problem 3

The left and right states of this Riemann problem are

\begin{matrix} (a_{g, L} = 0.55, u_{1, L} = 10.37 m / s, p_{L} = 80, 450 Pa) \\ (a_{g, R} = 0.55, u_{1, R} = 0.561 m / s, p_{R} = 24, 282 Pa) \end{matrix}

(41)

The solution consists of a left going shock wave, a contact wave, and a right going shock wave, as shown in Figure 4. The results obtained by different schemes at time $t = 1.0 s$ are shown in Figure 4. The results show that CE/SE scheme is more efficient than the modified CUP.

Figure 4.

Numerical results of Riemann problem by CE/SE and CUP schemes at t = 1.0 s.

Test problem 4

The left and right states of this Riemann problem are

\begin{matrix} (a_{g, L} = 0.35, u_{1, L} = 1.868 m / s, p_{L} = 192, 170 Pa) \\ (a_{g, R} = 0.30, u_{1, R} = 14.47 m / s, p_{R} = 196, 690 Pa) \end{matrix}

(42)

The solution composes of a left-going rarefaction wave, a contact discontinuity wave, and a right-going rarefaction wave, as shown in Figure 5. The solution profiles of phase velocities, volume fraction, and pressure are obtained from CE/SE and modified CUP schemes at time $t = 1.0 s$ , as shown in Figure 5. Once again, the CE/SE numerical scheme efficiently resolves the strong discontinuities as compared to the modified CUP scheme.

Figure 5.

Numerical results of shock tube problem by CE/SE and CUP schemes at $t = 1.0 s$ .

The last two problems were also considered in Fjelde and Karlsen.⁴ These test problems are considered the hardest for the numerical schemes. In these test problems, the source term $F_{w}$ is also involved. In these test problems, the gas slip-condition (7) is used with $C_{o} = 1.0$ and $u_{d} = 0.0 m / s$ . Here, the length of horizontal pipe is taken $1000 m$ with the diameter $0.1 m$ . In these test problems, the solution profiles are computed on 200 grid cells.

Test problem 5

The first $750 m$ part of pipe is filled with stagnant liquid, and the remaining area filled with 1% gas and 99% liquid. Initially, the pressure is $10^{5} Pa$ in the horizontal pipe. A pressure pulse is generated by increasing the inlet liquid rate from 0 to $0.3 kg / s$ in $0.0025 s$ . In the pure liquid region, the speed of sound is considered $1000 m / s$ and in the 1% gas region, the speed of sound is considered $100 m / s$ . Now, pressure waves move from left to right and will develop a strong shock wave, as shown in Figure 6. Because of sudden changes in the speed of sound and density of mixture, a large part of the pulse will be reflected. As a result of reflection, there will be a rarefaction wave. The solution profiles of liquid velocity and pressure at different times (t = 0.4, 0.9, and 1.3 s) are given in Figure 6, and the reference solutions are obtained for time t = 1.3 s. Both schemes behave very well for this hardest test problem, but the CE/SE numerical scheme shows better results at $t = 1.3 s$ .

Figure 6.

Numerical results of pressure pulse and liquid velocity at different times by CE/SE and CUP schemes.

Test problem 6

In this test problem, the pipe is initially filled with stagnant liquid. By injecting liquid and gas at the inlet, the mass flow rate of liquid as well as gas increased to $3.0$ and $0.02 kg / s$ , respectively, in $10 s$ . The pressure $p = 10^{5} Pa$ is kept constant at the outlet boundary. The solution profiles of liquid and gas mass flow rates at the four different positions $x = 250, 500, 750, and 1000 m$ are given in Figure 7, and the reference solutions are obtained at position $x = 1000 m$ .

Figure 7.

Numerical results at different positions by CE/SE and CUP schemes.

The gas and liquid passing through the horizontal pipe experience a decreasing pressure due to the frictional forces. As a consequence, the gas will expand, resulting in increased gas mass flow rates and the movement of liquid with larger velocity in front of the gas. Figure 7 shows that the liquid mass flow rate rises until the gas increases and decreases rapidly after passing the gas, while the rate of gas mass increases immediately when the gas arrives. The gas and liquid mass flow rates form sharp peaks at $x = 1000 m$ and then drop quickly. After $800 s$ approximately, the flow has stabilized. Once again, the results obtained by suggested numerical schemes are comparable and closely match with the results available in Fjelde and Karlsen.⁴ The CE/SE numerical scheme captures the sharp peak efficiently as compared to the modified CUP scheme, as shown in Figure 7.

Conclusion

In this article, the CE/SE scheme was extended to obtain the numerical solutions of the considered drift-flux model. The suggested numerical technique was capable to capture the volume-fraction discontinues without producing the oscillations. For comparison and validation, the modified CUP was also applied to solve the same drift-flux model. The number of test problems was considered. The numerical solutions obtained by CE/SE scheme verified the robustness, accuracy, and high resolutions for sharp discontinuities. A good agreement was found between the numerical solutions of both types of techniques, but CE/SE scheme has captured the strong shock waves and contact discontinuities more accurately and efficiently.

Footnotes

Handling Editor: James Baldwin

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

Rana Danish Aslam

References

Pelanti

LeVeque

. 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 Geo Res 1993; 98: 4231–4259.

Dufek

Bergantz

. Dynamics and deposits generated by the Kos Plateau Tuff eruption: controls of basal particle loss on pyroclastic flow transport. Geochem Geophy Geosy 2007; 8: 1–18.

Fjelde

Karlsen

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

Zuber

Findlay

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

Romate

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

Evje

Fjelde

. Hybrid flux-splitting schemes for a two-phase flow model. J Compu Phy 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 Flu 2006; 52: 679–705.

10.

Ishii

. Thermo-fluid dynamics of two-phase flow. Paris: Eyrolles.

11.

Evje

Fjelde

. Relaxation schemes for the calculation of two-phase flow in pipes. Mathema Comp Mod 2002; 36: 535–567.

12.

Evje

Fjelde

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

13.

Steinar

Tore

Munkejord

. A WIMF scheme for the drift-flux two-phase flow model, 2006, https://pdfs.semanticscholar.org/7738/53a3da1866491f69abe83b949087f4465cc3.pdf?_ga=2.13019440.648436184.1576746035-1540706140.1559042995

14.

Paula

Valdes

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

15.

Chang

. The method of space-time conservation element and solution element a new approach for solving the Navier-Stokes and Euler equations. J Compu Phy 1995; 119: 295–324.

16.

Chang

Wang

Chow

. New developments in the method of space-time conservation element and solution element-applications to two-dimensional time-marching problems (NASA-TM-106758), 1994, https://ntrs.nasa.gov/search.jsp?R=19950010470

17.

Chang

Wang

Chow

. The space-time conservation element and solution element method: a new high-resolution and genuinely multidimensional paradigm for solving conservation laws. J Compu Phys 1999; 156: 89–136.

18.

Liu

Wang

. The direct aero-acoustics simulation of flow around a square cylinder using the CE/SE scheme. J Algo Compu Technol 2009; 1: 525–537.

19.

Chang

Wang

. Application of the space-time conservation element and solution element method to one-dimensional convection-diffusion problem. J Compu Phys 2000; 165: 189–215.

20.

Loh

Hultgren

Chang

. Wave computation in in-compressible flow using the space-time conservation element and solution element method. AIAA J 2001; 39: 794–801.

21.

Loh

Hultgren

Chang

, et al. Noise computation of a shock-containing supersonic axisymmetric jet by the CE/SE method. In: 38th AIAA aerospace sciences meeting, Reno, NV, 10–13 January 2000, AIAA paper 2000-0475, https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20000025236.pdf

22.

Loh

Zaman

. Numerical investigation of transonic resonance with a convergent-divergent nozzle. AIAA J 2002; 40: 393–401

23.

Qamar

Mudasser

. On the application of a variant CE/SE method for solving two-dimensional ideal MHD equations. App Numer Math 2010; 60: 587–606.

24.

Wang

Chen

Liu

. The space-time CE/SE method for solving Maxwell’s equations in time-domain. In: 2002 IEEE international symposium on antennas and propagation national radio science meeting, San Antonio, TX, 16–21 June 2002. New York: IEEE.

25.

Shen

Wen

C-Y

Zhang

D-L

. A characteristic space-time conservation element and solution element method for conservation laws. J Compu Phy 2015; 288: 101–118.

26.

Shen

Wen

C-Y

. A characteristic space-time conservation element and solution element method for conservation laws II. Multidimensional extension. J Compu Phy 2016; 305: 775–792.

27.

Shen

Wen

C-Y

Liu

K-X

, et al. Robust high-order space-time conservative schemes for solving conservation laws on hybrid meshes. J Compu Phy 2015; 281: 375–402.

28.

Shen

Wen

C-Y

Parsani

, et al. Maximum-principle-satisfying space-time conservation element and solution element scheme applied to compressible multi-fluids. J Compu Phy 2017; 330: 668–692.

29.

Guan

Liu

Wen

C-Y

, et al. Numerical study on liquid droplet internal flow under shock impact. AIAA J 2018; 56: 3382–3387.

30.

Wen

C-Y

Massimi

Shen

. Extension of CE/SE method to non-equilibrium dissociating flows. J Compu Phy 2018; 356: 240–260.

31.

Fan

Guan

Wen

C-Y

, et al. Numerical study on the jet formation of simple-geometry heavy gas inhomogeneities. Phy Flu 2019; 31: 026103–026113.

32.

Kurganov

Lin

. On the reduction of numerical dissipation in central-upwind schemes. Commu Comp Phy 2007; 2: 141–163.

33.

Pico

Valdés

Ratkovich

, et al. Analysis of the drift flux in two-phase gas-liquid slug-flow along horizontal and inclined pipelines through experimental (non)-Newtonian and CFD Newtonian approaches. J Flu Flo Heat Mas Trans 2018; 5: 100–117.

34.

Kuila

Raja Sekhar

Zeidan

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