# Transient Simulation of Two-Phase Flows in Pipes

Q. H. Tran, D. Ferre, C. Pauchon, J. P. Masella

► To cite this version:

Q. H. Tran, D. Ferre, C. Pauchon, J. P. Masella. Transient Simulation of Two-Phase Flows in Pipes. Revue de l'Institut Français du Pétrole, 1998, 53 (6), pp.801-811. 10.2516/ogst:1998070. hal-02079048

HAL Id: hal-02079048

https://ifp.hal.science/hal-02079048

Submitted on 25 Mar 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
TRANSIENT SIMULATION
OF TWO-PHASE FLOWS IN PIPES

J.-M. MASELLA

Principia Recherche et Développement¹

Q. H. TRAN, D. FERRE and C. PAUCHON

Institut français du pétrole²

SIMULATION TRANSITOIRE DES ÉCOULEMENTS DIPHASIQUES DANS LES CONDUITES

La simulation transitoire des écoulements diphasiques gaz-liquide dans les conduites requiert des efforts de calcul considérables. Jusqu'à récemment, la plupart des codes de simulation commercialement disponibles étaient basés sur des modèles à deux fluides mettant en œuvre une équation de conservation de la quantité de mouvement pour chaque phase. Toutefois, dans les processus normaux de transport pétrolier dans les conduites, et spécialement de transport de pétrole ou de gaz, la réponse transitoire du système se révèle être relativement lente. Ainsi, il est raisonnable de penser qu'une forme simplifiée des équations de transport pourrait être suffisante pour reproduire de tels phénomènes transitoires. De plus, ce type de modèle pourrait être résolu numériquement au moyen d'algorithmes moins coûteux en temps calcul.

TRANSIENT SIMULATION OF TWO-PHASE FLOWS
IN PIPES

Transient simulation of two-phase gas-liquid flow in pipes requires considerable computational efforts. Until recently, most available commercial codes are based on two-fluid models which include one momentum conservation equation for each phase. However, in normal pipe flow, especially in oil and gas transport, the transient response of the system proves to be relatively slow. Thus, it is reasonable to think that simpler forms of the transport equations might suffice to represent transient phenomena. Furthermore, these types of models may be solved using less time-consuming numerical algorithms.

SIMULACIÓN TRANSITORIA
DE LOS FLUJOS DIFÁSICOS EN LOS CONDUCTOS.

La simulación transitoria de los flujos difásicos de gas líquido en los conductos, precisa esfuerzos de cálculo considerables. Hasta hace poco tiempo, la mayor parte de los códigos de simulación comercialmente disponibles se fundaban en modelos de dos fluidos que aplican una ecuación de conservación de la cantidad de movimiento para cada fase. No obstante, en los procesos normales de transporte de hidrocarburos por medio de conductos, y especialmente al tratarse del transporte de petróleo y de gas, la respuesta transitoria del sistema se manifiesta ser de forma relativamente lenta. Por todo ello, parece razonable pensar que una forma simplificada de las ecuaciones de transporte, podría ser de forma suficiente para la reproducción de tales fenómenos transitorios. Además, este tipo de modelo se podría resolver por vía digital mediante algoritmos menos costosos en cuanto a tiempos de cálculo.

<footnote>(1) ZI Brégaillon,
83857 La Seyne-sur-Mer</footnote>

<footnote>(2) 1 et 4, avenue de Bois-Préau,
92852 Rueil-Malmaison Cedex - France</footnote>
# INTRODUCTION

Transient simulations of gas-liquid flow in pipelines involve elaborate computer codes, the design and use of which demand tremendous effort. Several codes have been proposed so far, and supported by the oil and gas industry. Historically, OLGA is the first of these. Developed in Norway (Bendiksen et al. 1987, 1991), OLGA is a two-fluid model with an additional momentum equation for the droplet field. The PLAC code (Black et al. 1990), developed in England a few years later, is also a two-fluid model. At the Institut français du pétrole, another approach based on drift flux type models has been used. The resulting code TACITE is meant to become a commercial product (Pauchon et al. 1993, 1994). Another drift flux model, TRAFLOW, developed by the Shell Oil Company, should also be mentioned.

In the case of OLGA and PLAC, the model development and solution algorithms had been initiated earlier by the nuclear industry. In the nuclear industry, fast transients associated to Loss Of Coolant Accidents (LOCA) are of major interest, while in the oil and gas industry, the interest often lies in relatively slow transients, associated with the transport and subsequent release of slugs at receiving facilities. Under these conditions, one may consider the momentum equation to be a steady state force balance, thus leading to simpler and less elaborate calculations. Taitel et al. (1989) put forward a simplified model in which the gas mass flow rate is declared to be in steady state. Unfortunately, this approach lacked the ability to account for the transport time along the line for the case of gas flow rate variation at the inlet.

In this paper, we wish to study the behavior of various types of models under a few transient scenarios. Three different models were implemented, namely:

- A Two Fluid Model (TFM), based on one momentum conservation equation for each phase.

- A Drift Flux Model (DFM), based on one momentum conservation equation and an algebraic slip relation.

- A No Pressure Wave (NPW) model, based on an algebraic relation for the pressure drop and an algebraic slip relation.

It is to be expected that these three models have very different analytical properties. For instance, while TFM and DFM are hyperbolic models, the NPW model is a mixed parabolic/hyperbolic. Therefore, the numerical schemes in use will be different from one model to another.

First, we seek to analyze the differences in response due solely to the model equations. For the purpose of our study, the pipeline is visualized as a 1D element of length L. The coordinate along the pipe is called x. We also assume that the pipe properties such as inclination θ with respect to the horizontal, diameter D, roughness, etc. are constant with respect to x. Temperature is constant as well, and no mass transfer occurs between the two phases. The responses of the three transport models are compared to selected transient scenarios, which exemplify typical operational transient scenarios.

Then, we examine the three codes from the standpoint of computing efficiency. The latter reflects the trade-off between computing time and accuracy of the model. As for accuracy, it is defined in terms of the operating variables which are most significant for the end user, that is, the peak in outlet liquid flowrate subsequent to an increase in the inlet gas flowrate.

# 1 THE TWO-FLUID MODEL

In dispersed flow, the contrast between the two phasic velocities is small. Understandably, it is anticipated that the added value of the two-fluid model will be more apparent in the stratified regime. This is why we shall emphasize the stratified configuration in our presentation of TFM.

## 1.1 Transport Equations

TFM is governed by a set of four partial differential equations, the first two of which express mass conservation:

$$ \frac{\partial}{\partial t}[\rho_G R_G] + \frac{\partial}{\partial x}[\rho_G R_G V_G] = 0 \quad (1) $$

$$ \frac{\partial}{\partial t}[\rho_L R_L] + \frac{\partial}{\partial x}[\rho_L R_L V_L] = 0 \quad (2) $$

and the last two of which represent momentum balance:

$$ \frac{\partial}{\partial t}[\rho_G R_G V_G] + \frac{\partial}{\partial x}[\rho_G R_G V_G^2 + R_G \Delta P_G] + R_G \frac{\partial}{\partial x} P = \tau_G + \tau_i - \rho_G R_G g \sin \theta \quad (3) $$

$$ \frac{\partial}{\partial t}[\rho_L R_L V_L] + \frac{\partial}{\partial x}[\rho_L R_L V_L^2 + R_L \Delta P_L] + R_L \frac{\partial}{\partial x} P = \tau_L - \tau_i - \rho_L R_L g \sin \theta. \quad (4) $$
In Equations (3) and (4), P denotes the interface pressure and depends on masses ($\rho_G R_G$, $\rho_L R_L$), while $V_k$, $\rho_k$ and $R_k$ are respectively the velocity, the density and the volume fraction of phase $k \in \{G, L\}$. The variables $\tau_i$ and $\tau_k$ are the interfacial and wall momentum exchange terms, $\theta$ is the flowline inclination to the horizontal. Moreover we assume the constitutive law $R_G + R_L = 1$.

The quantities $\Delta P_G$ and $\Delta P_L$ correspond to the static head around the interface (De Henau and Raithby, 1995) defined as:

$$ \Delta P_G = P_G - P = -\rho_G \left[ \frac{1}{2} \cos\left(\frac{\omega}{2}\right) + \frac{1}{3\pi R_G} \sin^3\left(\frac{\omega}{2}\right) \right] gD \cos\theta \quad (5) $$

$$ \Delta P_L = P_L - P = \rho_L \left[ -\frac{1}{2} \cos\left(\frac{\omega}{2}\right) + \frac{1}{3\pi R_L} \sin^3\left(\frac{\omega}{2}\right) \right] gD \cos\theta \quad (6) $$

$\omega$ being the wetted angle.

<figure><img src="image_2.png" /><figcaption>Scheme 1

Definition of wetted angle.</figcaption></figure>

The densities are assumed to be given by $\rho_G = \frac{P}{a_G^2}$ and $\rho_L = \rho_L^0 + \frac{P - P_L^0}{a_L^2}$ where $a_k$ is the speed of sound in phase $k$, and $(\rho_L^0, P_L^0)$ is a reference state of the liquid phase.

It is convenient to rewrite Equations (1)-(4) in a more abstract way as:

$$ \frac{\partial}{\partial t}[\mathbf{w}] + \frac{\partial}{\partial x}[\mathbf{f}(\mathbf{w})] + \mathbf{r}(\mathbf{w})\frac{\partial}{\partial x}[P] = \mathbf{q}(\mathbf{w}) \quad (7) $$

where the vector:

$$ \mathbf{w}^T = (\rho_G R_G; \rho_L R_L; \rho_G R_G V_G; \rho_L R_L V_L) \quad (8) $$

represents conservative variables, as a function of which the vectors:

$$ \mathbf{f}^T = (\rho_G R_G V_G; \rho_L R_L V_L; \rho_G R_G V_G^2 + R_G \Delta P_G; \rho_L R_L V_L^2 + R_L \Delta P_L) \quad (9) $$

$$ \mathbf{r}^T = (0; 0; R_G; R_L) \quad (10) $$

$$ \begin{aligned} \mathbf{q}^T = & (0; 0; \tau_G + \tau_i - \rho_G R_G g \sin\theta; \\ & \tau_L - \tau_i - \rho_L R_L g \sin\theta) \end{aligned} \quad (11) $$

are computed. It is also possible to transform (7) into the quasi-linear form:

$$ \frac{\partial}{\partial t}[\mathbf{w}] + \mathbf{A}(\mathbf{w})\frac{\partial}{\partial x}[\mathbf{w}] = \mathbf{q}(\mathbf{w}) \quad (12) $$

Under appropriate conditions (Masella, 1997), e.g. if

$$ |V_G - V_L| < c_m \quad (13) $$

where $c_m$ is a pseudo-sound-velocity of the mixture, the matrix $\mathbf{A}(\mathbf{w})$ has four real eigenvalues:

$$ \lambda_1(\mathbf{w}) \le \lambda_2(\mathbf{w}) \le \lambda_3(\mathbf{w}) \le \lambda_4(\mathbf{w}) $$

as well as a base of eigenvectors. This property is commonly known as hyperbolicity. This hyperbolicity is depends on the relation (13). Note that, should the static head terms $\Delta P_G$ and $\Delta P_L$ be missing in Equations (3) and (4), then hyperbolicity would be lost as soon as $V_G \neq V_L$.

From the standpoint of physics, the extreme eigenvalues $\lambda_1$ and $\lambda_4$ are associated with acoustic waves and, therefore, can be very large, especially when the mixture is mainly composed of liquid. The intermediate eigenvalues $\lambda_2$ and $\lambda_3$ are associated with void fraction waves, and their orders of magnitude, which are close, are about those of fluid velocities.

## 1.2 Numerical Scheme

The resolution algorithm is based on a finite volume method, adapted to the non-conservative form (12). The reader is referred to Masella's Thesis (1997) for a more detailed description of the numerical implementation.

Let us divide the pipeline into a sequence of uniform cells $M_i = [x_{i-1/2}, x_{i+1/2}]$, the length of each of which is $\Delta x$. The unknowns are located at the center of the cells. Let $\mathbf{w}_i$ be the one associated to the cell $M_i$.
Discretization in space is based on a finite volume formulation of non-linear system (12) and reads

$$ \frac{d}{dt}[\mathbf{w}_i] + \frac{\mathbf{h}_{i+1/2} - \mathbf{h}_{i-1/2}}{\Delta x} + \mathbf{r}_i \frac{P_{i+1/2} - P_{i-1/2}}{\Delta x} = \mathbf{q}_i \quad (14) $$

where $\mathbf{h}_{i+1/2}$ and $P_{i+1/2}$ are obtained via a linearized Riemann problem at the interface $x_{i+1/2}$. This simplified Riemann problem, defined from the surrounding states $\mathbf{w}_i$ and $\mathbf{w}_{i+1}$, involves a rough Godunov solver and is referred to as the VFRoe solver (Masella et al., 1996).

Let $\mathbf{w}_{i+1/2}$ be the solution at the interface of the rough Riemann problem. Then $\mathbf{h}_{i+1/2}$ and $P_{i+1/2}$ are defined by:

$$ \mathbf{h}_{i+1/2} = \mathbf{f}(\mathbf{w}_{i+1/2}) \quad \text{and} \quad P_{i+1/2} = P(\mathbf{w}_{i+1/2}) \quad (15) $$

The quantity $\mathbf{r}_i$ can be put equal to $\mathbf{r}(\mathbf{w}_i)$.

We can note that the discretization of the non-conservative system (12) is obtained without any non-conservative reformulation like in Sainsaulieu 1995 and Toumi et al., 1996, but by a direct integration of the non-conservative system on each cell. A similar approach was used by Louis, 1995.

Note that the VFRoe scheme has been tested (Gallouet et al., 1996; Masella et al., 1996) with success on some hyperbolic conservative systems, such as Euler's and has given good results. These are quite comparable with the well-known Roe scheme.

The time discretization is explicit with respect to the intermediate eigenvalues and linearly implicit with respect to the extreme eigenvalues.

One of the trickiest problems here is to deal with the boundary conditions. Usually, two inlet mass flow rates and one outlet pressure are imposed. However, if $\lambda_1(\mathbf{w}) < 0 < \lambda_2(\mathbf{w}) \le \lambda_3(\mathbf{w}) \le \lambda_4(\mathbf{w})$ at the inlet, then a further inlet datum must be supplied. At first sight, a condition on the gas volumetric fraction $R_G$ or on the gas mass fraction $\alpha_G = \rho_G R_G / (\rho_G R_G + \rho_L R_L)$ seems to be relevant. Nevertheless, it comes as a surprise that such a choice does not always give rise to acceptable results, and in reality, the source terms must be taken into account.

As a matter of fact, most two-fluid type codes impose a von Neumann type boundary condition (De Henau and Raithby, 1995) such as:

$$ \frac{\partial}{\partial x}[R_L] = 0 \quad (16) $$

Others merely extrapolate the value of the liquid fraction. The staggered mesh approach used in OLGA,

with velocities defined at the mesh boundaries and pressure and void fraction defined within the cells, alleviates the need to specify an extra boundary condition.

A careful investigation into different terms in Equations (3) and (4) reveals that, in slow transients, the prevailing terms are the pressure drop and the source terms. When all remaining terms are neglected, Equations (3) and (4) degenerate into:

$$ \tau_G + \tau_L - (\rho_G R_G + \rho_L R_L) g \sin\theta - \frac{\partial}{\partial x}[P] = 0 \quad (17) $$

$$ R_G(\tau_L - \tau_i - \rho_L R_L g \sin\theta) \quad (18) $$

$$ -R_L(\tau_G + \tau_i - \rho_G R_G g \sin\theta) $$

This suggests to directly impose (18) as the third boundary condition at the inlet.

The next two figures illustrate what has just been said about boundary conditions. Two numerical runs were performed for a pipe of length $L = 10 000$ m. At the inlet, the flowrates are 10 kg/s (liquid) and 0.1 kg/s (gas). At the outlet, the pressure is maintained at 1 bar. Since three eigenvalues are positive at the inlet, a third boundary condition is required. We attempt to impose $R_G = 0.6$.

In the first run, the diameter is equal to $D = 0.25$ m, which corresponds to large friction terms. It can be seen from Figure 1 that the inlet value $R_G = 0.6$ is immediately dissipated and as a consequence, does not exercise any influence over the inner state of the mixture. This phenomenon is insensitive to the mesh size.

<figure><img src="image_3.png" /><figcaption>Figure 1</figcaption></figure>

Profile in space of the gas volumetric fraction for several mesh sizes (25, 50, 100 and 200 cells), with $D = 0.25$ m (zoom up to 5000 m).
In the second run, the diameter is increased to $D = 0.75$ m, so as to reduce the magnitude of friction terms. From Figure 2, it is now obvious that the inlet value $R_G = 0.6$ is perfectly compatible with the inner state.

<figure><img src="image_4.png" /><figcaption>Figure 2
Profile in space of the gas volumetric fraction for several mesh sizes (25, 50, 100 and 200 cells), with $D = 0.75$ m.</figcaption></figure>

## 2 THE DRIFT-FLUX MODEL

The Drift Flux Model is derived from the Two-Fluid Model by neglecting the static head terms $\Delta P_G$ and $\Delta P_L$ in the last two Equations (3) and (4) and replacing the two momentum equations by their sum. The main advantages of this 3-equation model are:

- The equations are in conservative form, which makes it easier to discretize by finite volume methods;

- The interfacial friction term $\tau_i$ is cancelled out in the momentum equation, although it appears in an additional algebraic relation called *slip law*;

- One does not have to work out a third boundary condition at the inlet;

- It is generally hyperbolic depending on the form of the *slip law*.

### 2.1 Transport Equations

Adding (3) and (4) together yields:

$$ \frac{\partial}{\partial t}[\rho_G R_G V_G + \rho_L R_L V_L] + \frac{\partial}{\partial x}[\rho_G R_G V_G^2 + \rho_L R_L V_L^2 + P] = \tau_G + \tau_L - (\rho_G R_G + \rho_L R_L) \sin \theta \quad (19) $$

The interfacial exchange term $\tau_i$ is no longer present in the above equation. This leads us to consider DFM, a new model which consists of three partial differential Equations, i.e. (1), (2) and (19). In order for DFM to be cast into the strictly conservative form:

$$ \frac{\partial}{\partial t}[\mathbf{w}] + \frac{\partial}{\partial x}[\mathbf{f}(\mathbf{w})] = \mathbf{q}(\mathbf{w}) \quad (20) $$

where the flux and the source:

$$ \mathbf{f}^T = (\rho_G R_G V_G; \rho_L R_L V_L; \rho_G R_G V_G^2 + \rho_L R_L V_L^2 + P) \quad (21) $$

$$ \mathbf{q}^T = (0; 0; \tau_G + \tau_L - (\rho_G R_G + \rho_L R_L)g \sin \theta) \quad (22) $$

are to depend only on:

$$ \mathbf{w}^T = (\rho_G R_G; \rho_L R_L; \rho_G R_G V_G + \rho_L R_L V_L) \quad (23) $$

it is necessary to introduce an algebraic relation called a *closure law* or *slip model*.

In stratified flow, the slip relation is obtained by combining the two momentum conservation Equations (3) and (4) in such a fashion that the pressure gradient vanishes. By neglecting derivatives with respect to time, static head and inertia terms in this combination, we end up with:

$$ \frac{\tau_G}{R_G} - \frac{\tau_L}{R_L} + \frac{\tau_i}{R_G R_L} + (\rho_G - \rho_L)g \sin \theta = 0 \quad (24) $$

For other flow regimes, the closure law may be much more sophisticated.

It is convenient to rewrite (19) in the quasi-linear form:

$$ \frac{\partial}{\partial t}[\mathbf{w}] + \mathbf{A}(\mathbf{w}) \frac{\partial}{\partial x}[\mathbf{w}] = \mathbf{q}(\mathbf{w}) \quad (25) $$

in which $\mathbf{A}$ is the Jacobian matrix of $\mathbf{f}$ with respect to $\mathbf{w}$. If hyperbolicity holds, there are three eigenvalues $\lambda_1 \le \lambda_2 \le \lambda_3$ and a base of eigenvectors. As is the case for TFM, the extreme eigenvalues are associated with acoustic waves, whereas the intermediate eigenvalue is associated with a void fraction wave.

### 2.2 Numerical Scheme

In this section, the basic ideas of the numerical scheme are outlined. Further details are to be found in Faille and Heintzé (1996).

Once again, let us divide the pipeline into regular cells $M_i = [x_{i-1/2}, x_{i+1/2}]$, whose size is denoted by $\Delta x$.
The unknowns are located at the center of the cells. Let $\mathbf{w}_i$ be the one associated with the cell $M_i$. Discretization in space reads:

$$ \frac{d}{dt}[\mathbf{w}_i] + \frac{\mathbf{h}_{i+1/2} - \mathbf{h}_{i-1/2}}{\Delta x} = \mathbf{q}_i \quad (26) $$

with:

$$ \mathbf{h}_{i+1/2} = \frac{1}{2}[\mathbf{f}(\mathbf{w}_i) + \mathbf{f}(\mathbf{w}_{i+1})] - \frac{1}{2}\mathbf{D}_{i+1/2}(\mathbf{w}_{i+1} - \mathbf{w}_i) \quad (27) $$

In this formula, $\mathbf{D}_{i+1/2}$ is a diffusion matrix that has to be built up from $\mathbf{w}_i$ and $\mathbf{w}_{i+1}$ to ensure sufficient stability. Extension to second order accuracy is achieved via the MUSCL strategy. Scheme (26)-(27) may appear to be more stable than standard Godunov-type schemes using approximate Riemann solvers (Godlewski and Raviart, 1996) at sides $x_{i+1/2}$ namely when some discontinuity, such as pipe slope change, at side $x_{i+1/2}$ occurs.

The time discretization is explicit with respect to the intermediate eigenvalue and linearly implicit with respect to the extreme eigenvalues.

# 3 THE NO-PRESSURE-WAVE MODEL

Experience with numerical simulations has shown that, in most transients of interest to the oil and gas transport industry, pressure waves do not have a strong effect on the initiation and transport of void waves. Hence, in a bolder step toward simplification, we would like to rule out the very existence of acoustic waves from the model equations.

## 3.1 Transport Equations

Preliminary numerical investigations into the order of magnitude of different terms in (18) suggest that, in slow transients, it is legitimate to neglect the inertia terms. Doing so amounts to replacing the mixture momentum Equation (19) by a local static force balance.

$$ \frac{\partial}{\partial x}[P] = \tau_G + \tau_L - (\rho_G R_G + \rho_L R_L)g \sin \theta \quad (28) $$

The NPW model is made up of (1), (2) and (28). The three partial differential equations are complemented by an algebraic slip law.

Let us introduce superficial velocities

$$ U_G = R_G V_G \quad (29) $$

$$ U_L = R_L V_L \quad (30) $$

$$ U_s = U_G + U_L \quad (31) $$

and let us consider

$$ \mathbf{v}^T = (R_G P, U_s) \quad (32) $$

as principal unknowns. This choice of variables turns out to be very handy for tackling the limiting case of incompressible flows. The equations of NPW can be summarized as:

$$ \frac{\partial}{\partial t}[\mathbf{e}(\mathbf{v})] + \frac{\partial}{\partial x}[\mathbf{f}(\mathbf{v})] = \mathbf{q}(\mathbf{v}) \quad (33) $$

with:

$$ \mathbf{e}^T = (\rho_G R_G; \rho_L R_L; 0) \quad (34) $$

$$ \mathbf{f}^T = (\rho_G R_G V_G; \rho_L R_L V_L; P) \quad (35) $$

$$ \mathbf{q}^T = (0; 0; \tau_G + \tau_L - (\rho_G R_G + \rho_L R_L)g \sin \theta) \quad (36) $$

The slip law is formally expressed as:

$$ U_G = \Psi(R_G, P, U_s) \quad (37) $$

where the function $\Psi$ may admit the pipe angle of inclination $\theta$ as a parameter.

Viviand (1996) shows that this model is a good approximation of the DFM as long as the phasic velocities are small compared to the sound wave velocities, which is true for most applications. He also proves that NPW always has a single finite eigenvalue, equal to:

$$ \lambda(\mathbf{v}) = \frac{\partial}{\partial R_G}[\Psi(\mathbf{v})] \quad (38) $$

The characteristic equation associated to this eigenvalue can be written as:

$$ \mathbf{I}^T(\mathbf{v}) \cdot \left\{ \frac{\partial}{\partial t}[\mathbf{v}] + \lambda \frac{\partial}{\partial x}[\mathbf{v}] \right\} = \mathbf{I}^T(\mathbf{v}) \cdot \mathbf{q}(\mathbf{v}) \quad (39) $$

$\mathbf{I}^T(\mathbf{v})$ being an appropriate left eigenvector. Besides, there exists an algebraically-double eigenvalue equal to $\infty$. On this ground, the model is qualified as mixed hyperbolic/parabolic. The number of characteristic equations associated with $\infty$ is usually one, i.e. (28), but may reach two for special thermodynamic laws.
## 3.2 Numerical Scheme

In this section, the basic ideas of the numerical scheme are outlined. Further details are to be found in Patault and Tran (1996).

The pipeline is divided into regular cells $M_i = [x_{i-1/2}, x_{i+1/2}$, the size of which is denoted by $\Delta x$. The unknowns are located at the center of the cells. Let $\mathbf{v}_i$ be the one associated to the cell $M_i$. Unlike the two previous schemes, the NPW equations are discretized in a centered way at the sides $x_{i+1/2}$. In other words:

$$ \frac{1}{2} \frac{d}{dt} [\mathbf{e}(\mathbf{v}_i) + \mathbf{e}(\mathbf{v}_{i+1})] + \frac{\mathbf{f}(\mathbf{v}_{i+1}) - \mathbf{f}(\mathbf{v}_i)}{\Delta x} = \frac{1}{2} [\mathbf{q}(\mathbf{v}_i) + \mathbf{q}(\mathbf{v}_{i+1})] \quad (40) $$

Time discretization of (40) is totally implicit. The final scheme is first order accurate in time, and second order in space.

Equation (40) has to be modified over a cell across which the sign of the eigenvalue $\lambda$ changes, to track shocks and rarefaction waves. An algorithm, inspired from Chattot and Mallet (1987), is set up to build the system of equations to be solved at each time step, according to the characteristic velocity sign configuration.

# 4 COMPARISONS OF THE MODELS

To compare the response of the three models, we first define a reference test case. This is done in Table 1. For this test case, we display the theoretical response of each transport model. By theoretical response, we mean the limit response obtained by decreasing the mesh size $\Delta x$ to 0.

<caption>TABLE 1</caption>

Definition of the reference test case

<table><tr><td>Pipe geometry</td><td>Length<br>Diameter<br>Configuration</td><td>420 m<br>3''<br>Horizontal</td></tr><tr><td>Transient scenario</td><td>Inlet gas flowrate<br>Inlet liquid flowrate<br>Outlet pressure</td><td>From 0.0048<br>to 0.0354 kg/s in 20 s<br>0.191 kg/s<br>1.68 bar</td></tr><tr><td>Fluid definition<br>and properties</td><td>Composition<br>Gas density<br>Liquid density<br>Gas viscosity<br>Liquid viscosity<br>Superficial tension</td><td>Air/Kerosen<br>2.418 kg/m³<br>813 kg/m³<br>0.76 10<sup>-5</sup> m²/s<br>0.22 10<sup>-5</sup> m²/s<br>0.03 N/m²</td></tr></table>

Next, we compare the response of the three codes to the experimental measurements. We also plot the computing time as a function of the accuracy for each model response.

Finally, we look at the response of NPW and DFM in a scenario giving rise to the severe slugging phenomenon. In this case, the differences in response are explained by arguments involving the transport models and the numerical schemes.

## 4.1 Transient Response of Models for the Reference Case

Let us start by describing the actual values of parameters for the reference case.

Figure 3 compares the transient response of DFM and NPW for $\Delta x = 2$ m. Below this value of the mesh size, no significant improvement in the accuracy of the solution for each model is observed. Therefore, we assimilate the curves to the “analytical” solutions of the models.

<figure><img src="image_5.png" /><figcaption>Figure 3
Comparison between DFM and NPW on outlet liquid flow rate.</figcaption></figure>

Note that the initial state is obtained solving each system TFM, DFM and NPW neglecting time differential terms and keeping boundary conditions constant.

To avoid discrepancies due to the slip model, we use in both cases a Zuber and Findlay (1965) slip model, which is an alternative to the slip model defined in Equation (24). The curves in Figure 3 exhibit very little difference in the transient response of the two models. This corroborates the theoretical prediction by Viviand, according to which pressure waves play a minor role in the propagation of the liquid.
<figure><img src="image_6.png" /><figcaption>Figure 4a
Comparison between DFM and TFM for outlet liquid flow rate</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4b
Comparison between DFM and TFM for liquid volume fraction.</figcaption></figure>

Figure 4 compares the transient responses of DFM and TFM, using a stratified slip model as specified in Equation (24). One must bear in mind the fact that TFM includes a pressure difference term due to the liquid fraction gradient, which is not included in DFM.

Although the transient imposed on the system is quite rapid (a 20-second ramp in inlet gas flow rate), the responses of DFM and TFM are quite similar. This result tends to indicate that TFM's decoupling of the phase velocities does not significantly affect the transient response of the model. In other words, it is the source terms that actually control the transient response.

<figure><img src="image_8.png" /><figcaption>Figure 5
Comparison of computing efficiency between DFM and NPW model.</figcaption></figure>

In Figure 5 we try to characterize the computing efficiency of the different codes. For this purpose, we define the “computing efficiency” as the relation between the computing time and the peak height in the liquid outlet flowrate. This maximum in the liquid outlet flowrate is important for pipeline design because it characterizes the outlet liquid flowrate surge into the separator following an increase in the inlet gas flowrate. The estimation of this liquid surge will eventually control the outlet separator volume design recommendation.

The diagram in Figure 5 shows that for a given accuracy, the mixed explicit/implicit solution used in the DFM code is faster than the implicit resolution in the NPW code.

## 4.2 Comparison with Experimental Data

Next we compare the code responses with experimental data (Vigneron, 1995) taken for the same test case. Figure 6 compares the experimental results with the transient response of the DFM using a slip model which takes into account the different flow regimes which can be encountered.

Figure 6a: Comparison between measured and computed liquid fraction using DF model.

The large amplitude oscillations observed during the transient period signal the occurrence of the slug flow regime. The recognition of this flow regime has important implications for the pressure drop calculation.
<figure><img src="image_9.png" /><figcaption>Figure 6a</figcaption></figure>

Comparison between measured and computed liquid fraction using DF model.

<figure><img src="image_10.png" /><figcaption>Figure 6b</figcaption></figure>

Comparison between measured and computed pressure using the DF model.

<figure><img src="image_11.png" /><figcaption>Figure 7</figcaption></figure>

Comparison between experimental and computed pressure using the DFM and TFM models with an imposed stratified configuration.

Figure 7 compares the same experimental results with the transient response of the DFM and TFM with the slip model assuming a stratified flow configuration.

From Figures 6b and 7, one can see the important effect of the flow regime prediction and subsequent slip relation on the pressure transient response of the model.

## 4.3 Severe Slugging Case

Lastly, we compare the DFM and the NPW models for a severe slugging test case.

Table 2 defines the severe slugging test case. Figure 8 shows the comparison between the two codes.

<caption>TABLE 2</caption>

Definition of the severe slugging test case

<table><tr><td>Pipe geometry</td><td>Length<br>Diameter<br>Configuration</td><td>60 m and 14 m<br>2''<br>Horizontal and vertical</td></tr><tr><td>Transient scenario</td><td>Inlet gas flowrate<br>Inlet liquid flowrate<br>Outlet pressure</td><td>0.000196 kg/s<br>0.07854 kg/s<br>1.0 bar</td></tr><tr><td>Fluid definition<br>and properties</td><td>Composition<br>Gas density<br>Liquid density<br>Gas viscosity<br>Liquid viscosity<br>Superficial tension</td><td>Air/Kerosen<br>1.0 kg/m³<br>1000 kg/m³<br>1.5 10<sup>-5</sup> m²/s<br>1.5 10<sup>-6</sup> m²/s<br>0.07 N/m²</td></tr></table>

Slight differences are observed in the period and amplitude of oscillations. However these can be attributed to the differences in treatment of the outlet boundary conditions. As can be seen, the DFM solution algorithm at the outlet does not allow for liquid going back into the pipeline during the reflux period. This problem is resolved by introducing a fictitious cell at the outlet in which we assume single phase gas flow.

## CONCLUSION

Extensive numerical tests have been carried out in order to compare the analytical response of the three following models:

- A Two Fluid Model with one momentum equation for each phase;
<figure><img src="image_12.png" /><figcaption>Figure 8a

Comparison between the DFM and the NPW models for pressure at riser bottom for a severe slugging case.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 8b

Comparison between the DFM and the NPW models for outlet liquid flow rate for a severe slugging case.</figcaption></figure>

- A Drift Flux Model with only one momentum equation for the mixture;

- A No Pressure Wave model reducing the momentum equation to a force balance.

On one hand, different hydrodynamic slip laws have been incorporated into these models. The results clearly show that the choice of the former has an critical effect on the dynamic response of the latter.

On the other hand, comparison between the different transport models confirms the fact that the source terms in the momentum equation do dominate the transient response of the model. This holds true for all models, but becomes overwhelmingly apparent in the case of TFM.

The computing efficiency as defined from the end user's perspective, that is, the peak in outlet liquid flow-rate subsequent to an increase in the inlet gas flowrate, has been analyzed for DFM and NPW. For a given accuracy, the DFM explicit/implicit solution is faster than the implicit solution employed in the NPW code.

## ACKNOWLEDGMENTS

We gratefully acknowledge the contributions of H. Viviand for his analysis of the NPW model and of S. Patault for its numerical implementation. We also wish to thank J. Brill (TUFFP) for allowing us to use one of his experimental data sets as a reference test case. This work was performed as part of the TACITE project supported by the Institut français du pétrole, Elf Aquitaine, and Total Exploration Production.

## REFERENCES

Bendiksen K., Brandt I., Jacobsen K.A. and Pauchon C. (1987) Dynamic simulation of multiphase transportation systems. *Multiphase Technology and Consequences for Field Development Forum*, Stavanger, Norway.

Bendiksen K., Malnes D., Moe R. and Nuland S. (1991) The dynamic two-fluid model OLGA: Theory and application. *SPE Production Engineering*, 6, 171-180.

Black P.S., Daniels L.C., Hoyle N.C. and Jepson W.P. (1990) Studying transient multiphase flow using the Pipeline Analysis Code (PLAC). *Journal of Energy Resources Technology*, 112, 25-29.

Chattot J.J. and Malet S. (1987) A box-scheme for the Euler equations. In: *Non-Linear Hyperbolic Problems*, Eds. Carasso C., Raviart P.A. and Serre D., *Lecture Notes in Mathematics*, 1270, Springer-Verlag.

De Henau V. and Raithby G.D. (1995) A transient two-fluid model for the simulation of slug flow in pipelines, 1-Theory. *Int. J. Multiphase Flow*, 21, 335-349.

Gallouët T, Masella J.M. (1996) On a rough Godunov scheme, *CRAS, Paris*, I-323, 77-84.

Faille I. and Heintzé E. (1996) Rough finite volume schemes for modelling two-phase flow in a pipeline. In: *Proceedings of the CEA, EDF, INRIA Course: "Méthodes numériques pour les écoulements diphasiques"*, INRIA, Rocquencourt, France.

Godlewski E. and Raviart P. A. (1996) *Numerical Approximation of Hyperbolic Systems of Conservation Laws*, Springer-Verlag, New York.

Louis X. (1995) *Doctoral Thesis*, University of Pierre et Marie Curie, Paris, France.

Masella J.M. (1997) *Doctoral Thesis*, University of Pierre et Marie Curie, Paris, France.

Masella J.M, Faille I. and Gallouët T. (1996) *On a Rough Godunov Scheme*, Submitted.

Patault S. and Tran, Q.H. (1996) Modèle et schéma numérique du code TACITE-NPW. *Technical Report* 42415, Institut français du pétrole, France.
Pauchon C., Dhulesia H., Lopez D. and Fabre J. (1993) TACITE: A comprehensive mechanistic model for two-phase flow. In: *Proceedings of the 6th International Conference on Multiphase Production*, Cannes, France.

Pauchon C., Dhulesia H., Binh-Cirlot G. and Fabre J. (1994) TACITE: A transient tool for multiphase pipeline and well simulation. *SPE Annual Technical Conference*, 28545, 25-28, New Orleans, LA, USA.

Sainsaulieu L. (1995); *Thèse d'habilitation*, University of Pierre et Marie Curie, Paris, France.

Taitel, Y, Shoham O. and Brill J. P. (1989) Simplified transient solution and simulation of two-phase flow in pipelines. *Chem. Eng. Sci.*, 44, 1353–1359.

Viviand H. (1996) Modèle simplifié d'écoulement dans les conduites pétrolières. *Technical Report RSY*, 43 069 02, Principia, France.

Vigneron F. (1995) Analysis of imposed two-phase flow transients in horizontal pipelines, Part 1-Experimental results. *Research Report*, Tulsa University Fluid Flow Project, USA.

Toumi I., Kumbaro A. (1996) An approximate linearized Riemann solver for a two-fluid model. *J. Comput. Physics.*, 124, 286-300.

Zuber N. and Findlay J.A. (1965) Average volumetric concentration in two-phase flow systems. *J. Heat Transfer Series* C 87, 453-458.

Final manuscript received in October 1998