Numerical investigations into a model of partially incompressible two-phase flow in pipes

TL;DR

This paper develops a numerical method for simulating two-phase liquid and gas flow in pipes, leveraging the independence of subsystems to simplify computations and improve efficiency.

Contribution

It introduces a staggered variable approach for separate numerical treatment of gas and liquid flows in a coupled two-phase model.

Findings

The staggered variable method simplifies the numerical solution process.

Using Roe flux for both subsystems yields accurate results.

Comparison with Nessyahu--Tadmor scheme demonstrates effectiveness.

Abstract

We consider a model for flow of liquid and gas in a pipe. We assume that the gas is ideal and that the liquid is incompressible. Under this assumption the resulting equations, expressing conservation of mass and momentum, splits into two subsystems such that the gas flow is independent of the liquid flow, and the liquid flow is described by a conservation law parametrized by the mass fraction of gas. When solving these equations numerically, we propose to stagger the gas and liquid variables with respect to each other. The advantage of this is that in finite volume methods one can use numerical flux functions designed for 2x2 systems of hyperbolic conservation laws to solve both the gas flow and the liquid flow, rather than a much more complicated numerical flux for the whole 4x4 system. We test this using the Roe numerical flux for both subsystems, and compare the results with results produced by using the second-order Nessyahu--Tadmor scheme for the second subsystem.

Tables

Table 1. (a) Liquid mass.

$N$	rel. $L^1$ error (Roe)	rel. $L^1$ error (NT)
$16$	$2.93$	$3.52$
$32$	$1.81$	$2.24$
$64$	$1.09$	$1.31$
$128$	$0.65$	$0.72$
$256$	$0.37$	$0.39$

Table 2. (a) Liquid mass.

$N$	rel. $L^1$ error (Roe)	rel. $L^1$ error (NT)
$16$	$2.93$	$3.52$
$32$	$1.81$	$2.24$
$64$	$1.09$	$1.31$
$128$	$0.65$	$0.72$
$256$	$0.37$	$0.39$

Table 3. (b) Liquid velocity.

$N$	rel. $L^1$ error (Roe)	rel. $L^1$ error (NT)
$16$	$10.55$	$12.26$
$32$	$6.46$	$8.72$
$64$	$3.90$	$5.13$
$128$	$2.34$	$2.81$
$256$	$1.34$	$1.48$

Table 4. (a) Liquid mass.

$N$	rel. $L^1$ error (Roe)	rel. $L^1$ error (NT)
$16$	$3.67$	$7.53$
$32$	$2.29$	$4.68$
$64$	$1.57$	$2.77$
$128$	$1.10$	$1.59$
$256$	$0.80$	$0.89$

Table 5. (a) Liquid mass.

$N$	rel. $L^1$ error (Roe)	rel. $L^1$ error (NT)
$16$	$3.67$	$7.53$
$32$	$2.29$	$4.68$
$64$	$1.57$	$2.77$
$128$	$1.10$	$1.59$
$256$	$0.80$	$0.89$

Table 6. (b) Liquid velocity.

$N$	rel. $L^1$ error (Roe)	rel. $L^1$ error (NT)
$16$	$6.78$	$11.16$
$32$	$4.42$	$8.44$
$64$	$2.92$	$5.76$
$128$	$1.99$	$3.09$
$256$	$1.42$	$1.70$

Equations

$$\text{partialderivative}t({\rho }_K{\alpha }_K)+\text{partialderivative}x({\rho }_K{\alpha }_K{v}_K)$$

$$\text{partialderivative}t({\rho }_K{\alpha }_K{v}_K)+\text{partialderivative}x({\rho }_K{\alpha }_K{v}_K^2+{\alpha }_K{p}_K)+{\rho }_K{\alpha }_K{g}_y\text{partialderivative}xh$$

$$s_K={\tau }_K\frac{{\sigma }_K}{A}\pm {\tau }_I\frac{{\sigma }_I}{A}-{\rho }_K{\alpha }_K{g}_x,$$

$${\alpha }_G+{\alpha }_L=1.$$

$${\alpha }_G{p}_G=\frac{m_G}{C_G}.$$

$$p_L$$

$$=\frac{m_G}{{\alpha }_G{C}_G}+\frac{1}{2}(1-{\alpha }_L)m_G+\frac{1}{2}\frac{m_L}{{\rho }_L}{m}_L$$

$$=\frac{m_G}{(1-{\alpha }_L)C_G}+\frac{1}{2}\left(1-\frac{m_L}{{\rho }_L}\right)m_G+\frac{m_L^2}{2{\rho }_L}$$

$$=\frac{m_G}{(1-\frac{m_L}{{\rho }_L})C_G}+\frac{1}{2}\left(1-\frac{m_L}{{\rho }_L}\right)m_G+\frac{m_L^2}{2{\rho }_L}.$$

$$\{\begin{array}{l}\text{partialderivative}t{m}_G+\text{partialderivative}x(m_G{v}_G)=0,\\ \text{partialderivative}t(m_G{v}_G)+\text{partialderivative}x\left(m_G{v}_G^2+\frac{m_G}{C_G}\right)=0,\end{array}$$

$$\{\begin{array}{l}\text{partialderivative}t{m}_L+\text{partialderivative}x(m_L{v}_L)=0,\\ \text{partialderivative}t(m_L{v}_L)+\text{partialderivative}x\left(m_L{v}_L^2+\frac{m_L{m}_G}{({\rho }_L-m_L)C_G}+\frac{m_L{m}_G}{2{\rho }_L^2}({\rho }_L-m_L)+\frac{m_L^3}{2{\rho }_L^2}\right)=0.\end{array}$$

$$u_t+f(u{)}_x$$

$$w_t+g(u,w{)}_x$$

$$u(x_j,t^n)$$

$$w(x_{j+\frac{1}{2}},t^n)$$

$$\text{partialderivative}tu+\text{partialderivative}xf(u)=0$$

$$u=\left(\begin{array}{c}m\\ mv\end{array}\right)\text{and}f(u)=\left(\begin{array}{c}mv\\ m{v}^2+p(m)\end{array}\right).$$

$$\text{partialderivative}tu+A(u^L,u^R)\text{partialderivative}xu=0,u(x,0)=\{\begin{array}{ll}{u}^L& x<0,\\ u^R& x>0.\end{array}$$

$$z(u)=\left(\begin{array}{c}{z}_1\\ z_2\end{array}\right)=\left(\begin{array}{c}\sqrt{m}\\ \sqrt{m}v\end{array}\right).$$

$$u(z)=\left(\begin{array}{c}{z}_1^2\\ z_1{z}_2\end{array}\right)\text{and}f(z)=\left(\begin{array}{c}{z}_1{z}_2\\ z_2^2+p(z_1^2)\end{array}\right)$$

$$\text{partialderivative}uz(z)=\left(\begin{array}{cc}2z_1& 0\\ z_2& z_1\end{array}\right)\text{and}\text{partialderivative}fz(z)=\left(\begin{array}{cc}{z}_2& z_1\\ 2z_2{p}^{\prime }(z_1^2)& 2z_2\end{array}\right)$$

$$\gamma (\xi )=z^L+(z^R-z^L)\xi ,\xi \in [0,1],$$

$$u^R-u^L={\int }_0^1\frac{du(\gamma (\xi ))}{d\xi }d\xi ={\int }_0^1\frac{du(\gamma (\xi ))}{d\gamma }{z}^{\prime }(\xi )d\xi ={\int }_0^1\frac{du(\gamma (\xi ))}{d\gamma }d\xi (z^R-z^L)=:B(z^R-z^L)$$

$$f(u^R)-f(u^L)={\int }_0^1\frac{df(\gamma (\xi ))}{d\xi }d\xi ={\int }_0^1\frac{df(\gamma (\xi ))}{d\gamma }{\gamma }^{\prime }(\xi )d\xi ={\int }_0^1\frac{df(\gamma (\xi ))}{d\gamma }d\xi (z^R-z^L)=:C(z^R-z^L)$$

$$f(u^R)-f(u^L)=A(u^R-u^L).$$

$$C:={\int }_0^1\frac{df(z(\xi ))}{d\gamma }d\xi =\left(\begin{array}{cc}\overline{z_2}& \overline{z_1}\\ 2\overline{z_1{p}^{\prime }(z_1^2)}& 2\overline{z_2}\end{array}\right)$$

$$B:={\int }_0^1\frac{du(z(\xi ))}{d\gamma }d\xi =\left(\begin{array}{cc}2\overline{z_1}& 0\\ \overline{z_2}& \overline{z_1}\end{array}\right)$$

$$B^{-1}=\frac{1}{2{\overline{z_1}}^2}\left(\begin{array}{cc}\overline{z_1}& 0\\ -\overline{z_2}& 2\overline{z_1}\end{array}\right)$$

$$\overline{h(z_k)}={\int }_0^{1}h(z_k^L+(z_k^R-z_k^L)\xi )d\xi ,k=1,2$$

$$A=C{B}^{-1}=\frac{1}{2{\overline{z_1}}^2}\left(\begin{array}{cc}0& 2{\overline{z_1}}^2\\ 2\overline{z_1}\overline{z_1{p}^{\prime }(z_1^2)}-2{\overline{z_2}}^2& 4\overline{z_1}\overline{z_2}\end{array}\right).$$

$${\lambda }_{1,2}=\frac{\overline{z_2}}{\overline{z_1}}\mp \sqrt{\frac{\overline{z_1{p}^{\prime }(z_1^2)}}{\overline{z_1}}},$$

$$|A|=R\mathrm{diag}(|{\lambda }_1|,|{\lambda }_2|)R^{-1}=\frac{1}{{\lambda }_2-{\lambda }_1}\left(\begin{array}{cc}{\lambda }_2|{\lambda }_1|-{\lambda }_1|{\lambda }_2|& |{\lambda }_2|-|{\lambda }_1|\\ {\lambda }_2{\lambda }_1^2-{\lambda }_1{\lambda }_2^2& {\lambda }_2^2-{\lambda }_1^2\end{array}\right)$$

$$F^{\text{Roe}}(u^L,u^R)=\frac{1}{2}\left(f\left(u^L\right)+f\left(u^R\right)\right)-\frac{1}{2}\left|A\left(u^L,u^R\right)\right|\left(u^R-u^L\right)$$

$$u_j^{n+1}=u_j^n-\frac{\Delta t}{\Delta x}\left(F^{\text{Roe}}\left(u_j^n,u_{j+1}^n\right)-F^{\text{Roe}}\left(u_{j-1}^n,u_j^n\right)\right).$$

Full text

1 Introduction

We consider the four-equation two-fluid model for stratified pipe flow with zero mass transfer. This model can be derived by averaging of the conservation equations across the cross-sectional area (see for example [11, 6, 3, 2]). The resulting equations are commonly written

[TABLE]

where $\rho_{K}$ denotes the density, $\alpha_{K}$ denotes the volume fraction and $v_{K}$ the velocity of phase $K$ ($K$ is “gas” or “liquid”). Furthermore, $h$ is the height of the interface between the two fluids and the momentum sources are

[TABLE]

where $\tau_{K}$ and $\tau_{I}$ are the wall respectively interface shear stress, $\sigma_{K}$ and $\sigma_{I}$ are the wetted lengths (see Figure 1) of the $K$ phase and the interface respectively and $A$ the pipe cross-sectional area. The axial and transverse components of the gravitational acceleration are given by $g_{x}=g\sin\phi$ and $g_{y}=g\cos\phi$, where $\phi$ is the pipe inclination relative to the horizontal plane. Figure 1 illustrates the pipe geometry and the aforementioned quantities.

Note that the mass $m_{K}=\rho_{K}\alpha_{K}$ and momentum $m_{K}u_{K}$ are conserved properties. We assume that phase $K$ is occupied by either gas, $K=G$, or liquid, $K=L$, and that the two phases are segregated from each other (see Figure 1). Note that for two phase flow

[TABLE]

In the following we will make some simplifying assumptions. First we assume zero momentum sources and no pipe inclination. Next we assume that the gas is ideal so that $\rho_{G}=p_{G}C_{G}$ for some constant $C_{G}$. Then

[TABLE]

Lastly, assuming that the liquid is incompressible, with constant density $\rho_{L}$, and assuming hydrostatic balance we find

[TABLE]

Therefore the system of equations we consider is

[TABLE]

If we let ${u}=(m_{G},m_{G}v_{G})$ and ${w}=(m_{L},m_{L}v_{L})$ we see that this can be written as

[TABLE]

for some nonlinear functions $f$ and $g$. This means that the $4\times 4$ model partially decouples into two $2\times 2$ systems, where the second system is dependent on the first, but not vice versa. The situation is reminiscent of so-called triangular systems of (scalar) conservation laws considered in [1], in which the convergence of a finite volume scheme was proved for such a system, in the case that ${u}$ and ${w}$ are scalars. In our case $u$ and $w$ are vectors, but when approximating, we can stagger the discretizations of $u$ and $w$ so that the approximation of $u$ is continuous (constant) across each cell interface in the approximation of $w$. The advantage of this is that the numerical flux for $w$ can be a standard numerical flux for $2\times 2$ system of conservation laws. In this paper we investigate using the Roe flux as a numerical flux for $w$. While not investigated in the present work, higher order methods can easily be built from first order numerical fluxes by higher order reconstruction.

The rest of this paper is organized as follows. In Section 2 we detail the construction of the Roe method for the gas and the liquid phase, as well as the Nessahu–Tadmor second order scheme. In Section 3 we describe how to calculate the exact solution to some Riemann problems, and test our scheme on these.

2 Numerical Methods

In this section we will describe the numerical methods we use to find approximate solutions to the two-phase flow problem (1.2). For the gas phase we will use Roe’s method and for the liquid phase we will use either Roe’s method or the second-order nonstaggered Nessyahu–Tadmor scheme.

We discretize the domain $[a,b]$ using two grids that are staggered with respect to each other. Let $\Delta x=(b-a)/N$ and $\Delta t=\frac{T}{M+1}$ and approximate $u=(m_{G},m_{G}v_{G})$ and $w=(m_{L},m_{L}v_{L})$ as follows:

[TABLE]

where $x_{j}=a+j\Delta x$, $j=0,\ldots,N$ and $x_{j+\frac{1}{2}}=a+(j+\frac{1}{2})\Delta x$, $j=0,\ldots,N-1$, as well as $t^{n}=n\Delta t$, $n=0,\ldots,M$. Thus the approximation to $u(\cdot,t^{n})$ is a piecewise constant function which may have discontinuities at $x_{j+\frac{1}{2}}$, and the approximation to $w(x,t^{n})$ is a piecewise constant function which may have discontinuities at $x_{j}$. Figure 2 illustrates the staggered grid and the approximations to $u$ and $w$.

2.1 The Roe method

Roe’s Riemann solver is a method to approximate the solution of the Riemann problem for the hyperbolic system

[TABLE]

with

[TABLE]

Here $m$ can be either $m_{G}$ or $m_{L}$ and $v$ can be $v_{G}$ or $v_{L}$ respectively, depending on which phase we consider. Roe’s approximate Riemann solver [10] (cf. also [4] and [8]) yields the solution of the following linear problem

[TABLE]

The so-called Roe matrix $A(u^{L},u^{R})$ has to satisfy the following properties:

1. (1)

$A(u^{L},u^{R})$ is diagonizable with real eigenvalues, 2. (2)

$A(u^{L},u^{R})\to f^{\prime}(u)$ as $u_{L},u_{R}\to u$ and 3. (3)

$A(u^{L},u^{R})(u^{R}-u^{L})=f(u^{R})-f(u^{L})$.

The first property is required for hyperbolicity and the second for consistency with the original conservation law. The third property ensures that single shocks of the nonlinear system are also discontinuous solutions of the linear system. In order to find such a matrix $A=A(u^{L},u^{R})$ we define the parameter vector

[TABLE]

Then $u$ and $f(u)$ can be written as

[TABLE]

and therefore

[TABLE]

Considering the straight line segment

[TABLE]

where $z^{L}=z(u^{L})$ and $z^{R}=z(u^{R})$, we find

[TABLE]

and similarly

[TABLE]

Then, by setting $A=CB^{-1}$ we have

[TABLE]

Note that

[TABLE]

and

[TABLE]

and thus

[TABLE]

where

[TABLE]

for a function $h$. Therefore

[TABLE]

The matrix $A$ has the eigenvalues

[TABLE]

which are real provided that $m>0$ and $p^{\prime}(m)>0$, and eigenvectors $r_{1,2}=\begin{pmatrix}1\\ \lambda_{1,2}\end{pmatrix}$. From this Roe’s method can be defined using

[TABLE]

where $R=(r_{1},r_{2})$ is the matrix consisting of the right eigenvectors to define the numerical flux

[TABLE]

which can then be used in the following Roe scheme:

[TABLE]

2.1.1 The Roe method for the gas phase

The pressure in the gas phase is

[TABLE]

with derivative

[TABLE]

Thus the Roe matrix $A$ becomes

[TABLE]

where

[TABLE]

The eigenvalues of $A(u^{L},u^{R})$ are

[TABLE]

and hence

[TABLE]

2.1.2 The Roe method for the liquid phase

The pressure in the liquid phase is

[TABLE]

with

[TABLE]

In order to use the Roe scheme here, we will have to evaluate the integral

[TABLE]

in (2.1) numerically. This can be done up to machine precision with an appropriate quadrature since $P_{m_{L}}$ is a rational function. After calculating the eigenvalues of the resulting matrix $A(u^{L},u^{R})$ given by (2.2) we can then assemble the Roe matrix $A(u^{L},u^{R})$ and the corresponding flux $F^{\text{Roe}}(u^{L},u^{R})$ as in (2.3) respectively (2.4).

2.2 Nonstaggered second-order Nessyahu–Tadmor scheme for the

liquid phase

We will now discuss a second-order scheme for the liquid phase. We want to solve the conservation law

[TABLE]

with

[TABLE]

and $u$ as before. Starting from the (staggered) second-order Nessyahu–Tadmor scheme [9]

[TABLE]

where $\lambda=\frac{\Delta t}{\Delta x}$ and the discrete derivatives $g_{j+\frac{1}{2}}^{\prime}$ and $w_{j+\frac{1}{2}}^{\prime}$ are made precise below. We use the averaging procedure described in [7] to get a nonstaggered version. To this end we reconstruct a piecewise-linear interpolant through the staggered cell-averages at time $t^{n+1}$:

[TABLE]

where the staggered discrete derivative $w^{\prime}_{j}$ is given by

[TABLE]

with

[TABLE]

Here, the minmod limiter is defined as

[TABLE]

We then average these interpolants over the cell $(x_{j},x_{j+1})$ to obtain a nonstaggered scheme:

[TABLE]

Here, $(w^{\prime}_{j})$ and $(w^{\prime}_{j+\frac{1}{2}})$ are given by (2.7) and

[TABLE]

respectively, where $\Delta w_{j}^{n}=w_{j+\frac{1}{2}}^{n}-w_{j-\frac{1}{2}}^{n}$. Note that, in our application, the flux $g$ depends on both $w$ and $u$, i.e.,

[TABLE]

where the discrete derivatives of the flux are given by

[TABLE]

and $u_{j+\frac{1}{2}}^{n}$ denotes the solution of the Riemann problem

[TABLE]

discussed in the previous section evaluated at $x=x_{j+\frac{1}{2}}$, i.e.,

[TABLE]

where $u_{k}^{1}$ and $u_{k}^{2}$ denote the first and second component of $u_{k}$ respectively (see [8]). Figure 3 illustrates all the values of $(u_{j}^{n})$ and $(w_{j+\frac{1}{2}}^{n})$ that are used when calculating $w_{j+\frac{1}{2}}^{n+1}$.

Note that the same averaging procedure detailed above can also be applied for higher-order schemes, cf. [7].

3 Numerical experiments

In order to have available exact solutions to which we can compare our approximations, we consider Riemann problems, i.e., the initial value problem for (1.2) where the initial data consists of a single jump between two constant values, viz.,

[TABLE]

If we fix the left state, we can find Riemann problems and their solution by following solution curves (rarefaction curves and Hugoniot loci) in phase space, see [5]. In the following sections we will first detail the ingredients necessary to find Riemann problems and their entropy solutions, and then present some example test cases and finally compare our schemes to those test cases.

3.1 Ingredients to solve the Riemann problem

We will now describe how one can find solutions to certain Riemann problems in two cases. Define

[TABLE]

and introduce the variables $q_{L}=m_{L}v_{L}$ and $q_{G}=m_{G}v_{G}$, and the flux function

[TABLE]

so that (1.2) reads

[TABLE]

The Jabobian of $F$ reads

[TABLE]

where $P_{m_{L}}=\frac{\partial P}{\partial m_{L}}$ and $P_{m_{G}}=\frac{\partial P}{\partial m_{G}}$. This matrix has eigenvalues

[TABLE]

Observe that this system is not strictly hyperbolic since the eigenvalues of the $\mu$ families can coincide with the eigenvalues of the $\lambda$ families.

3.1.1 Shocks

In the following, we will let

[TABLE]

denote the jump in some quantity $a$. Given $m^{L},m^{R}$ and either $v^{L}$ or $v^{R}$ the Rankine-Hugoniot loci in either phase can be computed by solving

[TABLE]

for $v^{R}$ respectively $v^{L}$, i.e. solving

[TABLE]

Here, the shock speed is

[TABLE]

Note that all four variables change over shocks of the $\lambda$ families. More specifically, rearranging the terms in (3.3) for the liquid phase we find

[TABLE]

Hence, using (3.5) and (3.6) in (3.3) we get

[TABLE]

Here, $\llbracket P(\cdot,\cdot)\rrbracket=P(m_{G}^{R},m_{L}^{R})-P(m_{G}^{L},m_{L}^{L})$.

Thus given a shock in the gas phase separating $(m_{G}^{L},v_{G}^{L})$ and $(m_{G}^{R},v_{G}^{R})$ traveling with speed $s$, and values $m_{L}^{L}$ and $v_{L}^{L}$ we can find the right state $(m_{L}^{R},v_{L}^{R})$ by solving

[TABLE]

for $m_{L}^{R}$ and computing the corresponding $v_{L}^{R}$ with (3.6).

3.1.2 Rarefactions of the $\mu$ families

Since the eigenvectors of the Jacobian (3.1) associated with $\mu_{1,2}$ are

[TABLE]

Rarefaction curves in the liquid phase are solutions of

[TABLE]

or

[TABLE]

If we use $v_{L}=v_{L}(m_{L},q_{L})$ instead, we get

[TABLE]

From this we see that rarefaction curves through a point $(m_{L}^{L},q_{L}^{L})$ are given by

[TABLE]

3.1.3 Rarefactions of the $\lambda$ families

The eigenvectors associated with $\lambda_{1,2}$, after renormalizing such that $\nabla\lambda_{1,2}\cdot r_{1,2}^{\lambda}=1$, are

[TABLE]

where

[TABLE]

By solving

[TABLE]

with the initial condition

[TABLE]

we get a rarefaction wave in all four variables.

3.2 Finding suitable Riemann problems and their solutions

We will now formulate two algorithms to find a Riemann problem whose solution consists of a shock in the liquid phase, followed by a shock over which all four components change, followed by a second shock in the liquid phase, and a Riemann problem whose solution consists of a rarefaction wave in the liquid phase, followed by a rarefaction wave over which all four components change, followed by a second rarefaction wave in the liquid phase respectively.

Algorithm 1: All-shock Riemann solution

1. (1)

Choose parameters $\begin{pmatrix}m_{G}^{L}\\ v_{G}^{L}\end{pmatrix},m_{G}^{R}$ as well as $m_{L}^{L},m_{L}^{\prime}$ and $m_{L}^{R}$. 2. (2)

Find {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}v_{G}^{R}} by solving the Rankine-Hugoniot equation (3.4) in the liquid phase such that the shock separating the left state $\begin{pmatrix}m_{G}^{L}\\ v_{G}^{L}\end{pmatrix}$ and right state \begin{pmatrix}m_{G}^{R}\\ {\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}v_{G}^{R}}\end{pmatrix} (with speed $s$) is a Lax shock of the $\lambda_{1}$ family. 3. (3)

Find {\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}v_{L}^{L}} by solving the Rankine-Hugoniot equation (3.4) in the liquid phase such that the shock separating the left state \begin{pmatrix}m_{L}^{L}\\ {\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}v_{L}^{L}}\end{pmatrix} and right state $\begin{pmatrix}m_{L}^{\prime}\\ v_{L}^{\prime}\end{pmatrix}$ (with speed $s_{L}^{L}$ such that $s_{L}^{L}<s$) is a Lax shock of the $\mu_{1}$ family. Here, $\llbracket P\rrbracket=p(m_{G}^{L},m_{L}^{\prime})-p(m_{G}^{L},m_{L}^{L})$. 4. (4)

Solve H({\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}m_{L}^{\prime\prime}},m_{L}^{\prime},m_{G}^{L},m_{G}^{R},v_{L}^{\prime},s)=0 for {\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}m_{L}^{\prime\prime}} and find the corresponding $v_{L}^{\prime\prime}$ with equation (3.6) so that the liquid phase has a shock with the same speed as the one in the gas phase. 5. (5)

Check that $s<\mu_{2}(m_{L}^{\prime\prime},v_{L}^{\prime\prime},m_{G}^{R})$. If not, choose different values for the parameters and start again. 6. (6)

Find {\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}v_{L}^{R}} by solving the Rankine-Hugoniot equation (3.4) such that the shock separating the left state $\begin{pmatrix}m_{L}^{\prime\prime}\\ v_{L}^{\prime\prime}\end{pmatrix}$ and right state \begin{pmatrix}m_{L}^{R}\\ {\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}v_{L}^{R}}\end{pmatrix} (with speed $s_{L}^{R}>s$) is a Lax shock of the $\mu_{2}$ family. Here, $\llbracket P\rrbracket=p(m_{G}^{R},m_{L}^{R})-p(m_{G}^{R},m_{L}^{\prime\prime})$.

Following this procedure we arrive at the following constant states in an all-shock solution of the Riemann problem:

[TABLE]

It is straightforward to extend this algorithm to include two shock waves associated with $\lambda_{1}$ and $\lambda_{2}$ by repeating Algorithm 1 for a shock of the $\lambda_{2}$ family and paying close attention to the ordering of the wave speeds.

Algorithm 2: All-rarefaction Riemann solution:

1. (1)

Chose parameters $\begin{pmatrix}m_{G}^{L}\\ v_{G}^{L}\end{pmatrix},\begin{pmatrix}m_{L}^{\prime}\\ v_{L}^{\prime}\end{pmatrix}$ and $m_{L}^{L}$ as well as $m_{L}^{R}$. 2. (2)

Find $\begin{pmatrix}m_{G}^{R}\\ v_{G}^{R}\end{pmatrix}$ and $\begin{pmatrix}m_{L}^{\prime\prime}\\ v_{L}^{\prime\prime}\end{pmatrix}$ by solving (3.8) to get a rarefaction over which all four components change of the $\lambda_{1}$ family. 3. (3)

Find {\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}v_{L}^{L}} by solving (3.7) such that we get the rarefaction wave connecting the left state \begin{pmatrix}m_{L}^{L}\\ {\color[rgb]{0,0.6,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.6,0}v_{L}^{L}}\end{pmatrix} and right state $\begin{pmatrix}m_{L}^{\prime}\\ v_{L}^{\prime}\end{pmatrix}$ is of the $\mu_{1}$ family. 4. (4)

Find {\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}v_{L}^{R}} by solving (3.7) such that we get the rarefaction wave connecting the left state $\begin{pmatrix}m_{L}^{\prime\prime}\\ v_{L}^{\prime\prime}\end{pmatrix}$ and right state \begin{pmatrix}m_{L}^{R}\\ {\color[rgb]{0,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{0,.5,.5}v_{L}^{R}}\end{pmatrix} is of the $\mu_{2}$ family.

Following this procedure we arrive at the following constant states in an all-rarefaction solution of the Riemann problem:

[TABLE]

It is straightforward to extend this algorithm to include two rarefaction waves associated with $\lambda_{1}$ and $\lambda_{2}$ or combine it with the algorithm above to include a $\lambda_{2}$ shock.

Experiment 1: All-shock Riemann solution

By following Algorithm 1 we can create our first test case:

[TABLE]

Here we have used $p(m_{G})$ and $p(m_{G},m_{L})$ as defined in (2.5) respectively (2.6) with $C_{G}=1$ and $\rho_{L}=1$. The solution to this Riemann problem consists of a shock in the liquid phase followed by a shock over which all four components change followed by a second shock in the liquid phase. Specifically, we have

[TABLE]

so all shocks are Lax shocks.

Experiment 2: All-rarefaction Riemann solution

Following Algorithm 2 we arrive at the second test case:

[TABLE]

Again we used $p(m_{G})$ and $P(m_{G},m_{L})$ as defined in (2.5) respectively (2.6) with $C_{G}=1$ and $\rho_{L}=1$. To solve the ODE (3.8) numerically we used a very high-order Runge–Kutta method and evaluated the integrals (3.7) using appropriate quadrature rules. Since

[TABLE]

the rarefaction waves are ordered as described above.

3.3 Comparison with numerical solutions

Figures 4 and 5 show the exact and numerical solutions of Experiment 1 and 2 respectively with open boundaries. Here we divided the spatial domain $[-5,5]$ into $54$ respectively $55$ ($N=55$) cells and use the time discretization parameter $\Delta t=\frac{1}{150}$ to approximate the solution at $T=1$.

Tables 1 and 2 show the relative $L^{1}$ errors between each of the two numerical approximations and the exact solution of the liquid mass and liquid velocity in Experiment 1 and 2 respectively for several values of the spatial grid size. We observe that the Roe scheme as an approximate Riemann solver has an overall smaller error in all instances.

4 Conclusion

We have considered a model for two-phase flow in pipes where we assumed hydrostatic balance, that the gas is ideal, and that the liquid incompressible. These assumptions allow us to partially decouple the $4\times 4$ system. Therefore, we can numerically treat the mass and momentum equation of the gas phase independently of the liquid phase. Staggering the numerical grid for the two subsystems with respect to each other further simplifies the numerical treatment since the approximation of the gas phase is constant at the cell interfaces of the approximation of the liquid phase.

As a starting point, we proposed to use the Roe scheme to numerically solve the gas subsystem and either the Roe scheme or the nonstaggered, second-order Nessyahu–Tadmor scheme for the liquid subsystem. In order to compare our numerical methods to exact solutions we further provided two classes of Riemann problems admitting only shocks respectively rarefactions. These classes can readily be modified to contain more than one shock, more than one rarefaction, or both shocks and rarefactions.

In our numerical experiments we compared numerical solutions of the liquid phase generated by the Roe scheme to the second-order Nessyahu–Tadmor scheme for certain Riemann problems.

Acknowledgements

We would like to thank Trygve Karper and Gunnar Staff from Schlumberger for helpful discussions on the topic.