A novel pressure-free two-fluid model for one-dimensional incompressible multiphase flow

A novel pressure-free two-fluid model formulation is proposed for the simulation of one-dimensional incompressible multiphase flow in pipelines and channels. The model is obtained by simultaneously eliminating the volume constraint and the pressure from the widely used two-fluid model (TFM). The resulting `pressure-free two-fluid model' (PF-TFM) has a number of attractive features: (i) it features four evolution equations (without additional constraints) that can be solved very quickly with explicit time integration methods; (ii) it keeps the conservation properties of the original two-fluid model, and therefore the correct shock relations in case of discontinuities; (iii) its solutions satisfy the two TFM constraints exactly: the volume constraint and the volumetric flow constraint; (iv) it offers a convenient form to analytically analyse the equation system, since the constraint has been removed. A staggered-grid spatial discretization and an explicit Runge-Kutta time integration method are proposed, which keep the constraints exactly satisfied when numerically solving the PF-TFM. Furthermore, for the case of strongly imposed boundary conditions, a novel adapted Runge-Kutta formulation is proposed that keeps the volumetric flow exact in time while retaining high order accuracy. Several test cases confirm the theoretical properties and show the efficiency of the new pressure-free model.


Introduction
The two-fluid model (TFM) for one-dimensional multiphase flow is an important model to study, for example, the behaviour of oil and gas transport in long pipelines. In this article we study its incompressible simplification. Like in single-phase flow models, the incompressibility assumption necessitates the use of careful space-and time-discretization methods. In the incompressible TFM, the main difficulty arises from the presence of the volume constraint (the two phases should together fill the pipeline), its associated volumetric flow constraint (the mixture velocity field should be divergence-free), and the role of the pressure.
One approach for numerically solving the incompressible TFM is to eliminate the pressure from the four-equation system (+ 1 constraint) and to rewrite this system into a two-equation system (+ 2 constraints). This leads to the 'no-pressure wave' model or the 'fixed-flux' model (FFM) suggested by Watson [22], and used for example in [1,6,12,14]. Another two-equation model is the 'reversed density' model developed by Keyfitz et al. [10] and employed for example in [21]. In these models, the pressure is generally computed as a post-processing step. A general problem with these two-equation systems is that they are only valid for smooth solutions. This means that in the presence of shocks different jump Email address: B.Sanderse@cwi.nl, corresponding author (B. Sanderse) conditions are obtained than for the four-equation model [1] (an example of shock waves in the two-fluid model are roll waves [22]). Furthermore, the fixed-flux assumption often limits these studies to stationary boundary conditions.
Another approach is to stick to the four-equation formulation and keep the pressure, and to use a pressure-correction method. A pressure equation is then typically obtained by substituting the momentum equations in the combined mass conservation equation, either with employing the volume constraint equation [11], or without employing the volume constraint equation [8,12,20]. In our recent work [17] we proposed a constraint-consistent pressure equation, which was shown to be both a generalized Riemann invariant of the continuous equations, and a hidden constraint of the semi-discrete differential algebraic equation system. The proposed Runge-Kutta time integration method was shown to be high-order accurate while satisfying both the volume constraint and the volumetric flow constraint. However, this approach is computationally relatively costly, as one needs to solve a pressure Poisson equation at each stage of the Runge-Kutta method, with the additional difficulty that the Poisson matrix depends on the actual solution and therefore changes in time.
In this paper we propose a novel incompressible two-fluid model formulation in which the pressure and the volume constraint have been eliminated. The elimination process will be such that the resulting pressure-free model still satisfies both the volume constraint and the volumetric flow constraint. Compared to the four-equation approach, our novel formulation is much faster to evaluate since it does not require the solution of a Poisson equation. Compared to the two-equation approach, the new model keeps the correct shock solutions, and furthermore it does not rely on the assumption of the 'fixed-flux'.
The outline of this paper is as follows: in section 2 the incompressible two-fluid model equations are given, in section 3 the new pressure-free model is derived, and in sections 4 and 5 appropriate spatial and temporal discretizations are proposed. Results are shown in section 6.

Governing equations for incompressible two-phase flow in pipelines
The incompressible two-fluid model can be derived by considering the stratified flow of liquid and gas in a pipeline (for a recent discussion of the two-fluid model, see for example [12]). The main assumptions that we make are that the flow is one-dimensional, stratified, incompressible, and isothermal. The transverse pressure variation is introduced via level gradient terms. Surface tension is neglected. This leads to the following system of equations, which we call the 'original' TFM [17]: where supplemented with the volume constraint equation: In these equations the subscript denotes either gas ( ) or liquid ( ). The model features four evolution equations, one constraint equation, and five unknowns (U and ), which are a function of the independent variables (coordinate along the pipeline axis) and (time). The primitive variables , etc. can be expressed fully in terms of U , e.g.
denotes the density (assumed constant), the cross-sectional area of the pipe, and (also referred to as the hold-ups) the cross-sectional areas occupied by the gas or liquid, the pipe radius, ℎ the height of the liquid layer measured from the bottom of the pipe, the phase velocity, the pressure at the interface, the shear stress (with the wall or at the interface), the gravitational constant, the local upward inclination of the pipeline with respect to the horizontal, and the components of gravity are = cos and = sin . See figure 1. The level gradient terms are given by [16]: for a pipe (compressible or incompressible flow); for incompressible flow in a channel one has instead In a pipe geometry, and ℎ are related by a non-linear algebraic expression, see e.g. [17]; in a channel one has = ℎ (assuming unit depth). The wetted and interfacial perimeters , and can be expressed in terms of the liquid hold-up or the interface height ℎ (see [17]). The source terms and constitute friction, gravity, and driving pressure gradient, given by The friction models are described in [17]. The source terms and do not contain spatial or temporal derivatives of the unknowns of the system (U and ).
Initial and boundary conditions determine whether the equations form a well-posed initial boundary value problem. It is well-known that this is not always the case: the two-fluid model equations can become ill-posed depending on the velocity difference between the phases [3,13,18]. In this paper we will restrict ourselves to test cases for which the model is well-posed. Furthermore, as will become clear in the next section, we require the initial conditions to be consistent, meaning that at = 0 the solution should not only satisfy the volume constraint but also the volumetric flow constraint:

The hidden constraints of the incompressible TFM
In this section we show how the pressure and the constraint can be removed from the incompressible TFM. For this purpose, we use the characteristic analysis of the incompressible two-fluid model proposed in [17]. In that work, it was shown that the volume constraint, equation (3), could be rewritten in two alternative forms.
The first form is the volumetric flow constraint, obtained by premultiplying the TFM equations with leading to The last equation is the volumetric flow constraint. Note that in this derivation we used the time-derivative of the volume constraint, so the volumetric flow constraint can be interpreted as a combination of mass equations and volume constraint. Integrating the volumetric flow constraint in space gives The second form is the pressure equation, obtained by premultiplying the TFM equations with leading to At the heart of this derivation is the time-derivative of equation (15), which in itself was obtained from the time-derivative of the volume constraint, see equation (14). In other words, equation (17) can be interpreted as an equation for the second-order time-derivative of the constraint: Note that the pressure gradient equation (17) is different from the expression commonly used in literature (see e.g. [7]), in which the gas and liquid momentum equations are added without scaling by their respective densities. Although such a pressure equation can be useful for postprocessing computations, it cannot be used to replace the volume constraint, because the volume constraint has not been employed in its derivation. In contrast, our derived pressure equation (17) is a combination of the momentum equations, mass equations, and the volume constraint.

The pressure-free TFM
The key insight of this paper is the following. The expression for the pressure gradient, equation (17), can be substituted back into the TFM equations. For this purpose, first write the pressure equation in short notation as wherê The pressure gradient term as present in the TFM then reads with The important conclusion is that the pressure gradient term in the TFM can be written in terms of a linear combination of the flux vector and the source terms. Upon substituting the expression for h(U ) into the TFM, we obtain the new pressure-free TFM (PF-TFM) as

PF-TFM:
where Equation (23) is a closed system of four equations in four unknowns (provided thaṫ ( ) is known): both the constraint and the pressure have been removed. The volume constraint is 'built into' this system of equations (it was used in the derivation of the pressure equation), and does not need to be provided as a fifth equation. We note that making the model pressure-free comes at a small 'price': one needṡ ( ) as an input to the system of equations. In other words: removing the pressure requires the prescription oḟ ( ). Fortunately,̇ ( ) equals zero in many cases, e.g. in case of steady inflow conditions with prescribed liquid and gas flow rates. For unsteady inflow, e.g. an increasing liquid and gas production,̇ ( ) is known from the boundary condition specification. Only in special cases, such as periodic boundary conditions, ( ) is unknown; we will assume it to be zero in that case and investigate by comparing PF-TFM solutions to TFM solutions to what extent this is indeed true.

Properties of the PF-TFM
The PF-TFM shares the property of the TFM that it cannot be written in full conservative form. In the TFM this manifests itself in the presence of the hold-up fractions in front of the pressure gradient term; in the PF-TFM this enters through the presence of the A(U ) matrix. Note that A(U ) is a singular matrix, so the system cannot be simplified by pre-multiplying with the inverse of A. The eigenvalues of A(U ) f (U ) U are exactly the same as those obtained from the analysis of the original TFM in primitive variables in [17], see Appendix A, so the PF-TFM has the same characteristic velocities as the TFM.
The PF-TFM does not require the solution of a Poisson equation for the pressure, which significantly accelerates the code (this will be shown in the results section). Given a solution U to the PF-TFM, the pressure can be obtained as a post-processing step by solving equation (19). An existing TFM implementation (based on a pressure Poisson equation) can easily be adapted by simply scaling the computed fluxes f using the weights from the A matrix, and adding thė ( ) term. The pressure Poisson equation can then be skipped.
Compared to the fixed-flux model (FFM), there are two important advantages. First, in contrast to the FFM, the PF-TFM retains the correct shock relations (i.e., those of the TFM). Second, the FFM is typically used with 'fixed-flux' assumption ( ) = constant [12], whereas this assumption is not made in our model.
Although the constraint is 'built into' the pressure-free TFM, upon discretization care should be taken, because the volume constraint is not built into the PF-TFM in its original form, equation (3), but in a differentiated sense, namely equation (18). This will be discussed next.

Spatial discretization on a staggered grid to preserve constraints
In the derivation of the pressure-free model, equation (23), the following important property was silently used: the term present in the time-derivative in the momentum equations, e.g. 3

=
, is the same term as in the flux term in the mass equations. This property needs to be satisfied at a discrete level in order to obtain a pressure-free TFM that satisfies the volume constraint; this is achieved here by the use of a staggered grid.
The staggered grid consists of = 'pressure' and = + 1 'velocity' volumes (with slight changes depending on the boundary conditions). Although the pressure is not present in the new formulation, this concept is so well-known that we stick to the terminology. The midpoints of the velocity volumes lie on the faces of the pressure volumes. The pressure, hold-up and phase masses are defined in the centre of the pressure volumes, whereas the velocity and momentum are defined in the centre of the velocity volumes. For details we refer to [16]. The unknowns are the vector of conservative variables ( ) ∈ ℝ 2 +2 , which include the finite volume sizes (we use the notation to distinguish from the vector U that appears in the continuous equations): Including the volume size in the vector of unknowns leads to a clear physical interpretation (e.g. 1 has units of mass) and allows generalization to the case of time-dependent finite-volume sizes. 1 , 2 ∈ ℝ and 3 , 4 ∈ ℝ contain the values of mass and momentum at the mass and velocity volumes, respectively, and are only a function of time. With these vectors, the expression for the gas hold-up and gas velocity is (similar to (4)): where Ω ∈ ℝ × is a diagonal matrix that contains the pressure volume sizes, Ω ∈ ℝ × a diagonal matrix containing the velocity volume sizes, and ∈ ℝ × an interpolation matrix from pressure to velocity points.
,1∕2 ,3∕2 ,1 = 1 We start with conservation of mass for the gas phase. Integration of the first equation in (23) indirection over a pressure volume gives: A crucial detail in this equation is that the convective fluxes are directly expressed in terms of the momentum 3 . For all volumes, the discretization is written as where is a spatial differencing matrix, see [17]. For clarity: 1 denotes the first component of the spatial discretization of A f + AS + ċ .
For conservation of momentum we proceed in a similar way. Integration of the third equation in (23) in -direction over a velocity volume gives: where the notation indicates an entry in the A matrix. For the approximation of the convective and level gradient terms present in , see [17]. In summary, the discretization of the momentum equation of the gas phase is written as and similar for 4 . We stress that in contrast to the original TFM, 3 does not contain the pressure, but instead contains the A-matrix and volumetric source term involvinġ . The entire spatial discretization (mass and momentum equations) can be summarized as We check whether the volume and volumetric flow constraints are satisfied with this spatial discretization. The first-order time-derivative of the volume constraint is given by The second-order time-derivative then follows as The last term can be simplified by realizing that the A matrix in the PF-TFM has been constructed such that: wherė ∈ ℝ = 1̇ ( ) is a vector with ( ) in each entry. This expression can be easily checked by substituting the expressions for 3 and 4 (see equation (29)), and for A and c, see equation (24). See also equation (B.12). Consequently, equation (33) gives The last equality holds because application of the differencing matrix to a constant vector yields 0. Integrating equation (35) twice in time, given an initial condition that satisfies the volume constraint ( ( 0 ) + ( 0 ) = ), then leads to In other words, the volume constraint remains satisfied in time by employing a spatial discretization on a staggered grid.

Basic formulation
Whereas the previous section showed that the pressure-free TFM posed requirements on the spatial discretization in order to be constraint-consistent, no further requirements are necessary for the temporal discretization: a simple forward Euler time discretization already keeps the constraint property. The extension to generic explicit Runge-Kutta methods is then straightforward. Note that due to the absence of the pressure, the spatially discretized system can be advanced straightforwardly in time (without the need to solve a Poisson equation).
The forward Euler step is given by An important term present in ( ) is the discrete approximation to the last term in equation (29), being the volumetric flow source term. We split ( ) as The forward Euler discretization of the last term, as given by equation (37), yields: An alternative formulation (option B) of the source term will be discussed in section 5.2. As before, the key question is whether the constraint is satisfied at the new time level +1 (assuming that it was satisfied at the current time level ): Here we used option A for the discretization of the volumetric source term. Given an initial condition for 3 and 4 at = 0 which satisfies the volumetric flow constraint, and given that ̇ ( ) = 0, the volumetric flow constraint can be written as a telescopic sum over all previous time steps: The volume constraint at the new time level, equation (40), can thus be evaluated as In summary, when applying the forward Euler method to equation (31) and starting from consistent initial conditions, both the volume constraint and the volumetric flow constraint are satisfied at each time step. Note that neither equation (41), nor equation (43), guarantees that the solution +1 will satisfy the exact flow rate ( +1 ), since only the time-derivativė ( ) has been used. However, when instead of option A, we apply the following first-order discretization of the volumetric flow source term option B: into equation (41), not only the volumetric flow constraint is satisfied, but in addition the actual value of the volumetric flow stays equal to the specified volumetric flow when marching in time: In other words, when adapting the forward Euler discretization for the volumetric flow source term, it is possible to satisfy the volumetric flow constraint and the actual value of ( ) exactly. The generalization of this idea to high-order time integration methods requires careful considerations, as will be discussed next.

High-order accuracy and integration oḟ ( )
The high-order time integration methods we consider are explicit Runge-Kutta (RK) methods. Explicit RK methods form an excellent combination of accuracy and stability for convection-dominated systems and, in contrast to the standard TFM, can be applied to the proposed PF-TFM in a straightforward manner, because the constraint has been eliminated. When considering option A for handling the volumetric flow source term, the resulting method reads: where and are the coefficients of the Butcher tableau, and the number of stages. Since explicit RK methods can be seen as a combination of forward Euler steps, it is straightforward to prove that they satisfy the two constraints, both for the stage values , and for +1 . When considering option B for handling the volumetric source term, there is an open question at which time level to evaluate ( ) at the intermediate stages. One possibility that leads to an exact volumetric flow is to take at each stage ( )( ( , ) − ( )).
However, this is choice does not fit in the Runge-Kutta framework, and the resulting method suffers from loss of order of accuracy. The question how to employ option B in conjunction with a high-order method can be answered with the following insight. The time discretization corresponding to option B, equation (45), is obtained when deriving the pressure-free model on the fully discrete level instead of on the continuous level. In other words: option A is obtained by first eliminating the pressure from the two-fluid model, and then discretizing the resulting system; option B is obtained by first discretizing the two-fluid model, and then eliminating the pressure. The details of this derivation are shown in Appendix B. With this insight, a high-order RK method for option B is straightforward to construct: eliminate the pressure from the (fully discrete) high-order RK methods proposed in [17], and rewrite as a time integration method for the PF-TFM.
We focus on the case = 3, i.e. a three-stage third-order method, for which the following method results: It can be checked through substitution that the volumetric flow remains equal to the specified volumetric flow ( ), both at the stage levels and at the new time level. It is important to remark that, although this is strictly speaking not a Runge-Kutta method applied to (31), it is a Runge-Kutta method applied to the original TFM. The order of accuracy of this method was analyzed in [17], and a specific third-order method was proposed: This method is such that the lower diagonal entries ( 21 , 32 ) are nonzero (this is required in equations (52)-(53)), and order reduction that can result from unsteady boundary conditions is avoided (details in next section). The -coefficients are used to determine the intermediate time levels ( = + Δ ) and should not be confused with the vector c or as present in the PF-TFM.
In summary, we have proposed high-order time integration of the PF-TFM for both 'weak' (option A) and 'strong' (option B) imposition of the volumetric flow source term. In both cases the two constraints are satisfied, and in the latter case also the volumetric flow remains equal to the specified value. The second approach requires a specialized Runge-Kutta method, which is slightly more involved to implement. Note that in the special case thaṫ ( ) = 0, we havê = and option B becomes a classic Runge-Kutta method, equivalent to option A. In this case we will resort to a classic four-stage fourth order Runge-Kutta method.

Boundary conditions
The boundary condition treatment follows the characteristic approach outlined in [17], which has the advantage that it does not require the pressure in the formulation. Here we shortly summarize the (important) case of unsteady inflow conditions, with prescribed gas and liquid (mass) flows, ( ) and ( ), respectively. There are two choices for the boundary conditions: prescription in strong form (prescribe the actual value of and at each time step), or in weak form (prescribe the time-derivativeṡ anḋ and integrate in time with the Runge-Kutta method).
It is important to note that the boundary condition prescription should be consistent with the discretization of the volumetric source term ( )̇ ( ) in order to make sure that remains uniform in space, so that the volumetric flow constraint is satisfied. The weak boundary condition prescription is consistent with option A. The strong boundary condition prescription is consistent with option B.

Keeping the constraint exact in time by preventing machine error accumulation
In deriving equation (40) the error in the volume constraint at the previous time step was imposed to be exactly equal to zero. In actual computations, the value of + − might not exactly equal zero due to the presence of machine precision errors. By simply leaving this term in the numerical algorithm, one can avoid accumulation of machine precision errors that could spoil the accuracy of the computations [5,19]. In other words, one should not impose the volume constraint term in (40) to be zero, but instead simply compute its value based on the solution at the last time step ( + − ) and use this while time stepping. In practice we achieve this by computing the error in the volume constraint after each stage of the Runge-Kutta method: and adding this term to the solution after each stage: where Note that this does not affect 3 and 4 , so it does not change the volumetric flow.

Kelvin-Helmholtz instability
We perform the classic viscous Kelvin-Helmholtz instability test case, see e.g. [11,17] and table 1 for parameter values. We choose values for and , = 1 m∕s and = ∕ = 0.9 and then compute the gas velocity and the body force necessary to sustain the steady solution: The velocity difference − is below the 'Kelvin-Helmholtz' instability limit [11,17], which means that, at least initially, the initial boundary value problem is well-posed. In addition, the conditions are such that the solution is linearly unstable, so we can study the growth of waves. We impose thaṫ = 0, We perturb the steady state by imposing a sinusoidal disturbance with wavenumber = 2 and a small amplitude. Linear stability analysis [11,16] gives the following angular frequencies : The negative imaginary part of 2 indicates the instability in the solution. We choose the initial perturbation related to 1 to be zero, implying that a single wave with frequency 2 results, and the linearized solution is where ΔU is obtained by choosing the liquid hold-up fraction perturbation aŝ = 10 −3 , and then computing the perturbations in the gas and liquid velocity from the dispersion analysis [11]. The initial condition which follows by taking = 0 s does in general not satisfy (11) exactly, and we therefore perform a projection step to make the initial conditions consistent, see [17], in such a way that the volumetric flow rate stays exact. Figure 3 shows the solution at = 1.5 s computed with = 40 finite volumes and Δ = 1∕100 s. The amplitude has grown with respect to the initial condition due to the negative imaginary part in 2 . At the same time, non-linear effects are causing wave steepening, which will lead to shock formation at later times. Figure 4 shows the error in the solution upon time step refinement, where the error is computed with respect to a reference solution computed with Δ = 1 × 10 −4 s. Figure 5 shows that the errors related to the volume constraint, the volumetric flow constraint, and the actual volumetric flow are satisfied until machine precision. Note that these errors are defined at each time instant as: volume constraint error: volumetric flow constraint error: volumetric flow error: and the maximum is taken over all elements in the vector.

Discontinuous solutions: roll waves
In this second test case we simulate roll waves on a periodic domain in order to (i) investigate the difference in ( ) between the proposed PF-TFM and the original TFM for the case of periodic boundary conditions, (ii) show that the PF-TFM captures the same discontinuous solution as the TFM, and (iii) give an indication of the computational speed-up that can be achieved with the PF-TFM.
Roll waves can form when simulating the Kelvin-Helmholtz instability of section 6.1 further into the nonlinear regime. The test case that we perform here is inspired by the experimental work of Johnson [9] and simulations of Akselsen [2], Holmås [7], and Sanderse [15]. The parameter settings are reported in table 2. We first compute a steady state solution: we prescribe the flow rates ∕ = 3.5 m∕s, ∕ = 0.35 m∕s, which leads to = 0.190 and d 0 d = −155.919 Pa∕m. The initial condition follows by applying a sinusoidal perturbation of the form (61), with a hold-up amplitude of 0.01, = 2 ∕ and = 4.597 − 0.068 chosen such that a single, unstable, wave is triggered (similar to the Kelvin-Helmholtz test case in section 6.1). One important difference with respect to the previous test case is that the interfacial friction factor in the expression for the interfacial stress ( = 1 2 | − |( − )) is taken as = ⋅ , with = 12.5. Another important difference is that the discretization of the convective terms in the momentum equation is performed with a first-order upwind scheme in order to prevent numerical oscillations at the discontinuity (shock wave), instead of the central scheme used in the first test case. Figure 6a shows the hold-up profile at = 100 s computed with both the PF-TFM and the TFM. At this time instant the wave has travelled approximately 73 times through the domain, covering a distance of about 219 m. The seemingly large difference in the position of the discontinuity (about 0.4 m) between the PF-TFM prediction and the TFM prediction is therefore in fact less than 0.2% of the total travelled distance. This difference is consistent with the relative difference in the volumetric flow rate predicted by the PF-TFM and TFM, which is also approximately 0.2% (see figure 6b). This difference reaches a constant value after approximately 50 seconds, as the roll wave then reaches a stationary solution. Note that the TFM prediction, which does not require a prescribed value of ( ), should be interpreted as being the most accurate solution. The PF-TFM prediction can be seen as a highly accurate approximation to this solution at a much reduced computational cost.
For this particular case ( = 320), the computational cost of the PF-TFM prediction is about 40% lower than of the TFM. For other grid sizes a similar cost reduction is observed, as is shown in figure 7. Figure 7a shows a quadratic scaling of CPU time with , which is caused by the fact that the simulations are run at a fixed CFL number (Δ = 1∕ ). Figure 7b shows that the obtained reduction is relatively independent of the number of grid cells. However, we should note that other factors can have an important influence, such as the type of pressure solver used in the TFM (here a direct solver was used) or the type of friction model. For example, when a more computationally expensive friction model like the Biberg model [4] is used, the gain in computational cost obtained with the PF-TFM compared to the TFM is probably less pronounced, because the contribution of the pressure Poisson solve to the total computational cost is smaller. On the other hand, the computations are performed here on a single CPU. Much larger speed-ups with the PF-TFM can be expected in a parallel implementation: the PF-TFM framework is fully explicit and can be easily parallelized, whereas the original TFM requires has an implicit component (the solution of a Poisson equation) for which good parallel scaling will be more difficult to obtain.

Hold-up wave propagation
The final test case considers the propagation of a hold-up wave through a horizontal pipeline caused by an increasing inlet gas flow rate, inspired by the test case in [14]. The parameters of the problem are described in table 3. The initial conditions are steady state production with inlet mass flows of liquid and gas of = 1 kg∕s and ,start = 0.02 kg∕s. Furthermore, we specify the following time-dependent gas flow rate at the inlet, which is such thaṫ (0) = 0 and for > 0̇ ( ) (and higher order derivatives) are nonzero. The volumetric flows ( ) and are given by ( )∕ and ∕ , respectively. The solutions with the weak BC prescription (option A) and with the strong BC prescription (option B) at = 1000 s with = 40 and Δ = 10 s are indistinguishable, see figure 9a. Figure 9b shows that the errors in volume constraint and volumetric flow constraint (equations (62)-(63)) remain at machine precision for both methods. In addition, option B has the property that the numerically computed volumetric flow rate remains exactly equal to the specified value. In option A, a small error is present, which decreases with third order upon time step refinement.
A more quantitative assessment is obtained by computing a reference solution at = 1000 s with = 40, Δ = 1 × 10 −2 s, and weak boundary conditions (option A). At this small time step, the volumetric flow error stays at machine precision also with weak imposition of and boundary conditions. Next we compute solutions at a sequence of much larger time steps, Δ = 40, 20, 10, …, both with weak BC (option A) and strong BC (option B). Figure 8 shows that both options converge according to the design order of accuracy (third order).

Conclusions
In this article we have proposed a novel pressure-free incompressible two-fluid model for multiphase pipeline transport, in which both the pressure and the volume constraint have been eliminated. Together with a staggered-grid spatial discretization and an explicit Runge-Kutta time integration method, this yields a new highly efficient method for solving the incompressible two-fluid model equations. The explicit nature also allows for a straightforward parallelization. Furthermore, the absence of the pressure make the model more amenable to analysis. For example, the PF-TFM offers a new viewpoint in studying the issue of the conditional well-posedness of the TFM. The 'price' to pay compared to the conventional form of the incompressible two-fluid model, is that the time-derivative of the volumetric flow rate needs to be specified as an additional input to the equations. Fortunately, in most cases for pipeline transport design (e.g. steady production, production ramp-up) this value is known from the boundary conditions.
A natural question is whether the original TFM 'pressure Poisson' formulation [17] is still useful, given that a more efficient pressure-free alternative has been proposed in this work. The answer is yes, for two reasons. First, the original formulation also works wheṅ is unknown (e.g. in the case of periodic boundary conditions). Second, as indicated in [17], the original formulation has the potential to be extended to multi-dimensional problems, with application to for example the volume-of-fluid approach. This is probably not possible with the proposed approach: the fact that the pressure gradient can be expressed directly in terms of the conservative variables by using the constraint is presumably only possible when dealing with one-dimensional problems. Extension of the pressure-free model to one-dimensional incompressible three-phase flow is straightforward. On the other hand, extension to compressible flow is not straightforward: in compressible flow, differentiation of the volume constraint leads to a pressure equation that involves the time derivative of the pressure. Consequently, the pressure gradient term in the momentum equations cannot be eliminated in the same way as in the incompressible model. These eigenvalues are the same as those of the TFM, see e.g. [17].

Appendix B. Pressure-free model from fully discrete two-fluid model equations
In this appendix we derive a discrete PF-TFM based on a fully discrete approximation of the TFM. This corresponds to 'option B'. The forward Euler discretization of the TFM reads [17]: ∈ ℝ × is an interpolation matrix, ∈ ℝ × a differencing matrix (similar to ), and ⊙ denotes elementwise multiplication. To arrive at the pressure-free model based on the discrete two-fluid model, we follow exactly the same steps as in the continuous case, outlined in section 3. We substitute the mass equations in the constraint and apply the expression for TFM, 12 : Note that the first term in brackets (.) corresponds to the matrix element 33 and the second term in brackets to 34 .
In summary, the PF-TFM forward Euler discretization derived from the fully discrete forward Euler discretization of the TFM reads +1 = + Δ ̂ ( ) + ( )( − −1 ), (B. 13) and is such that the volumetric flow stays exact when integrating in time (option B). The extension to a generic explicit Runge-Kutta method is a straightforward substitution exercise.