Energy-consistent formulation of the pressure-free two-fluid model

The pressure-free two-fluid model (PFTFM) is a recent reformulation of the one-dimensional two-fluid model (TFM) for stratified incompressible flow in ducts (including pipes and channels), in which the pressure is eliminated through intricate use of the volume constraint. The disadvantage of the PFTFM was that the volumetric flow rate had to be specified as an input, even though it is an unknown quantity in case of periodic boundary conditions. In this work, we derive an expression for the volumetric flow rate that is based on the demand for energy (and momentum) conservation. This leads to PFTFM solutions that match those of the TFM, justifying the validity and necessity of the derived choice of volumetric flow rate. Furthermore, we extend an energy-conserving spatial discretization of the TFM, in the form of a finite volume scheme, to the PFTFM. We propose a discretization of the volumetric flow rate that yields discrete momentum and energy conservation. The discretization is extended with an energy-conserving discretization of the source terms related to gravity acting in the streamwise direction. Our numerical experiments confirm that the discrete energy is conserved for different problem settings, including sloshing in an inclined closed tank, and a traveling wave in a periodic domain. The PFTFM solutions and the volumetric flow rates match the TFM solutions, with reduced computation time, and with exact momentum and energy conservation.

minimized. Important applications include oil and gas or CO 2 transport through pipelines 2,3 and nuclear reactor safety analysis. 4 We consider the incompressible form of the TFM, for (separated) stratified flow, assuming a hydrostatic balance between the fluids. 5 One reason this model is of interest, is its potential to dynamically simulate the transition from stratified flow to slug flow. 6,7 This stands in contrast to for example the drift-flux model, where the two velocities are not modelled separately: instead a closure relation is introduced for their difference. 8,9 The transition to slug flow can occur through hydrodynamic instabilities that arise naturally for the TFM. [10][11][12] However, the physical instability of the model is connected to a loss of hyperbolicity, an issue for which different interpretations and circumvention strategies have been proposed. [13][14][15][16][17] This issue is shared with the two-layer shallow water equations (TLSWE), as described for example by Chiapolino and Saurel. 18 Both models include one mass conservation and one momentum conservation equation for each fluid (totaling four conservation equations). While the TLSWE describe open channel flow, the TFM describes flow in a duct that is closed at the top. This introduces a constraint to the model, named the volume constraint, describing the demand that the two fluids must together fill the duct. While in the TLSWE the pressure is completely determined by the ambient pressure and the hydrostatic pressure, in the closed duct of the TFM the pressure can vary independently as a fifth variable, and a Poisson equation for it follows from the volume constraint and the derived volumetric flow constraint, which describes the incompressible nature of the flow. Furthermore, the TFM describes flow in ducts with arbitrarily shaped cross-sections.
The constraints and the nondynamic pressure variable complicate the numerical solution of the TFM, and its (nonlinear) stability analysis. The four equations of the TFM can be combined into two to form the fixed flux model (FFM), by taking the sum of the mass equations, and the difference of the momentum equations, which eliminates the pressure. 12,[19][20][21] In the process, smoothness of the solution is assumed, and as a result the model has incorrect jump conditions for shocks. 22 In addition, the volumetric flow rate is often assumed to be constant in time, which excludes unsteady inflow boundary conditions, and which is incorrect in case of periodic boundary conditions. Similar two-equation models have been given by Keyfitz et al. 23 and Etrati et al., 24 and a comparable approach has been taken for a two-layer shallow water model with rigid-lid assumption, 25 which describes flow in a closed channel. A different two-equation two-fluid model was derived by Jones and Prosperetti, 26 in which the volumetric flow constraint is used to eliminate both the pressure and the mass equations. Like the FFM, this model assumes smoothness of the solution, and the volumetric flow rate needs to be prescribed.
In a recent article by Sanderse et al., 27 it was shown how a four-equation pressure-free two-fluid model (PFTFM) can be derived, by eliminating the pressure using the volume constraint, but otherwise leaving the model in its original conservative form. This model requires no constraints, exhibits the correct shock relations, and allows a time-dependent volumetric flow rate.
An open issue with the PFTFM was that the volumetric flow rate needed to be prescribed, although for periodic boundary conditions it is a priori not known what value to take: it is a free parameter. Without an appropriate prescription for this free parameter, the solutions of the TFM and the PFTFM, when starting from identical initial conditions, will diverge in time, 27 which means that the PFTFM solution will be incorrect. In this work, the demand for energy conservation is used to find an appropriate expression for the rate of change of the volumetric flow rate. Since it is known that the TFM solution conserves energy, 28 replicating this property in the PFTFM naturally yields a PFTFM solution consistent with that of the TFM.
Due to the manner in which the PFTFM is derived from the TFM, it is natural to apply mass-and momentum conserving discretizations of the TFM to the PFTFM (unlike for two-equation models which are not formulated as mass and momentum conservation equations per fluid). In the present paper, we will show that, starting from our recently presented staggered grid energy-conserving discretization of the TFM, 28 an energy-conserving semi-discrete form of the PFTFM can be derived. This includes an expression for the rate of change of volumetric flow rate, as a function of the current solution throughout the domain. This discretization can be compared to energy-conserving discretizations of the shallow water equations, on a staggered grid, 29 or on a collocated grid. 30,31 These discretizations improve the stability of numerical solutions.
We extend the discretization with a streamwise potential energy conserving discretization of the streamwise gravity source terms, to make the discretization more complete and extend its applicability to more general cases. The discretization of these terms is similar to an existing staggered-grid energy-conserving discretization of terms in the SWE that describe the influence of a nonflat bottom topography. 29 However, in the SWE, the bottom topography is added to a level domain, while in the TFM the spatial coordinate follows the inclination of the duct. This paper is set up as follows. Section 2 introduces the TFM and the PFTFM, and shows how the latter is derived from the former. It points out the step in the derivation where, if the PFTFM's free parameter is not set correctly, the models could diverge. In Section 3, we analyze the energy behavior of the continuous PFTFM, as a means to find an expression for the free parameter, which removes the discrepancy between the TFM and the PFTFM. Additionally, it is demonstrated that the streamwise gravity terms can be included in an energy-conserving framework. Section 4 is the semi-discrete counterpart of Section 2, and derives the semi-discrete PFTFM from the semi-discrete energy-conserving TFM, before discussing both models' incorporation of the flow constraints. Section 5 uses semi-discrete energy conservation as a guide to find a semi-discrete expression for the free parameter, and introduces the energy-conserving discretization of the streamwise gravity terms. Finally, in Section 6, we present the results of numerical experiments, which show that with our proposed expression for the free parameter, and with our proposed discretization, the PFTFM solution matches that of the TFM perfectly.

Introduction to the TFM
The one-dimensional incompressible two-fluid model for stratified flow in channels and pipes can be derived by setting up integral mass and momentum balances for two control volumes: one spanning a section of the duct occupied by the upper fluid, and the other spanning a section of the duct occupied by the lower fluid (see Figure 1). Taking the limit in which the length of the control volumes goes to zero yields a set of equations for upper and lower fluid cross-sections A U and A L , and the cross-sectionally averaged velocities u U and u L . 15,32 These are functions of time t and the axial coordinate s. In this derivation, a few assumptions are made; one of the most important being that variations along the streamwise direction typically occur over distances larger than the duct diameter, resulting in an assumption of hydrostatic balance in the normal direction. 19,33 The cross-sectionally averaged equations can be written in terms of the conservative variables q(s, t): 28,34 q with The cross-sections are multiplied by the upper and lower fluid densities, U and L , in order to obtain a mass per unit length, and the velocities are multiplied by the cross-sections and densities to obtain a momentum per unit length. In (1), the fluxes are given by HereĤ U =Ĥ U (q) andĤ L =Ĥ L (q) are geometric terms which are measures for the potential energy relative to the interface height, defined in Appendix A. The gravitational acceleration in the normal direction is represented by g n = g cos( ), with the angle at which the duct is inclined.
F I G U R E 1 A schematic of stratified two-fluid flow in ducts (a circular pipe segment is shown as an example) described by the The pressure p appearing in the equations is the pressure at the interface. It is preceded by the vector d, which is given by Neglecting friction and surface tension, the source terms are given by ] .
Gravity in the streamwise direction is represented by g s = g sin( ) (see Figure 1), and a driving pressure gradient p 0 ∕ s may be imposed in case of periodic boundary conditions. The system of equations (1) is subject to the volume constraint which through differentiation to t and substitution of the mass conservation equations can be converted into the volumetric flow constraint: with Since Q = Q(q) and q = q(s, t), Q depends on t through q (but not on s due to (4)). We can use these constraints to set up an equation for the pressure. Typically, this equation is obtained by summing the momentum equations as follows: 35 and substituting (5) and using (4) to get withQ Finally, taking the derivative of (6) to s and applying (4) gives a Poisson equation for the pressure, which can be used instead of the volume constraint to close the model, acting as the fifth equation of the TFM (in addition to (1)):

Pressure-free two-fluid model
The PFTFM was introduced by Sanderse et al. 27 It is a form of the incompressible two-fluid model from which the pressure has been eliminated. To this end, (7) is multiplied by d and rewritten as with Substitution of (10) in (1) then yields the PFTFM Here, the term kQ (appearing in (10)) has been replaced by kK. This has been done because in this model,K =K(t) is a free parameter that needs to be set. We cannot use the local definition in terms of q ′ ∕ t (given by (8)), since substituting this definition would result in a singular system. 1 Instead, the rate of change of the volumetric flow rate needs to be specified explicitly, as a function that is independent of s. This directly enforces the volumetric flow constraint (see Section 2.4), without need for the pressure and its Poisson equation. We have given q ′ an apostrophe in order to distinguish the solution of the PFTFM from the solution of the TFM, given by q. The transformation from TFM to PFTFM will normally not change the solution (meaning that the original solution q(s, t) will still satisfy the transformed equations), unlessK(t) is set such that it differs from theQ(t) corresponding to q(s, t). In this case, with both models starting from the same initial conditions q(s, t 0 ) = q ′ (s, t 0 ), the solutions to the TFM and PFTFM will diverge over time so that q(s, t) ≠ q ′ (s, t).
The derivation described in this subsection relies on the existence of an s-independent value for the volumetric flow rate, which is guaranteed by (4) for this one-dimensional model. In multidimensional incompressible flow, the divergence-free constraint does not determine a unique volumetric flow rate, so this method cannot be applied; see Lteif 36 for a similar conclusion regarding the TLSWE with a rigid lid.

Boundary conditions
The boundary conditions are important to consider due to their relevance to the volumetric flow rate. At the boundaries, q 3 and q 4 are set, and the remaining variables follow via characteristic analysis. 35 This holds both for the TFM and the PFTFM. Consider a domain bounded by s 1 and s 2 . For a closed domain, the boundary conditions are given by q 3 (s 1 , t) = 0, q 4 (s 1 , t) = 0, q 3 (s 2 , t) = 0, and q 4 (s 2 , t) = 0. This directly determines Q = Q(q) at the boundaries through (5), and through (4) this determines Q throughout the domain.
For inflow boundary conditions, q 3 and q 4 at s 1 and s 2 may be (known) functions of time. Therefore, Q may be a function of time (through q(s, t)), but it is again fully determined by the boundary conditions. Through (4), Q will be determined for the full domain, also at an outflow boundary where q 3 and q 4 depend on the solution in the interior.
In case of periodic boundaries, the boundary conditions are q 3 (s 1 , t) = q 3 (s 2 , t) and q 4 (s 1 , t) = q 4 (s 2 , t). While this illustrates that Q at s 1 will be equal to Q at s 2 , these boundary conditions do not determine the actual value of Q. As a result, Q can vary freely with time. The degree to which it varies depends on the specifics of the initial conditions and the parameter values. For the TFM, Q can be seen as a component of the solution, for which the model equations must be solved.

2.4
The problem of setting the volumetric flow rate For the PFTFM, it is necessary to set the (rate of change of the) volumetric flow rate as an input to the system of equations through the free parameterK =K(t). This is the cost of eliminating the pressure p = p(s, t) as a variable. 27 By taking the dot product of (12) with l and substituting we find thatQ which makes clear that the value we set for K determines the actual volumetric flow rate of the PFTFM solution. This holds throughout the domain, which directly implies that the PFTFM solution satisfies the volumetric flow constraint (4), if the initial conditions satisfy the volumetric flow constraint. Similarly, take the dot product of (12) with with A defined as a function of q ′ by (3). If care is taken that the initial conditions satisfy the constraints, this equation demonstrates that the volume constraint is also naturally satisfied by the PFTFM. Due to the discussion in Section 2.3, for closed and inflow boundaries it is clear how to setK, but for periodic boundaries it is a priori unclear. In related models, the problem of setting the volumetric flow rate for periodic domains is usually avoided by assumingQ ′ = 0. 12,19 The same assumption has been made (at a later stage in the derivation) by Sanderse et al., 27 where in case of periodic boundaries the settingK = 0 has been used. However, this simplification yields differences between PFTFM and TFM solutions: we haveK(t) =Q ′ (t) ≠Q(t) and q ′ (s, t) ≠ q(s, t). Considering the TFM to be the ground truth, this renders PFTFM solutions incorrect. In this work we will introduce an expression forK based on the demand for momentum and energy conservation, that resolves this discrepancy.

Energy conservation for the continuous model
As stated above, an expression for the free parameterK can be determined based on the demand for energy conservation.
To this end, we will first give an outline of the steps required to prove that the mechanical energy is a secondary conserved quantity of the PFTFM, based on our proof of energy conservation for the regular TFM. 28 In general, the process involves defining an expression for the mechanical energy, and (based on this expression) manipulating the governing equations of the model to obtain a local energy conservation equation, from which global energy conservation trivially follows. Given an expression for the mechanical energy e(q), we can define the vector of entropy variables, following Fjordholm et al., 37 as Taking the dot product of this energy with the governing equations given by (12) yields where we have omitted the apostrophes, but emphasize that q refers to the PFTFM solution. The first term of this equation which makes (14) an equation for e∕ t. Now, if the remaining terms can be written in conservative form, meaning that they can be written as the one-dimensional divergence of some yet to be defined energy fluxes h and j (it will be determined in the current and the following subsection that this is indeed possible for the proposed model): then (14) reduces to the local energy conservation equation showing how the energy e(s, t) at a specific point in space changes due to an inflow or outflow, or due to source terms. Therefore, (15) and (16) are the conditions which need to be verified in order to prove energy conservation for the PFTFM.
In the case of periodic boundaries the difference between an energy flux at left boundary and the energy flux at the right boundary must be zero, and in case of closed boundaries the boundary values of the energy fluxes must themselves be zero. Therefore, for both periodic and closed boundaries, integrating (17) over the domain yields the global energy conservation equation Since the PFTFM is directly derived from the original TFM and describes the same physics, we expect that it conserves the same energy as the original TFM. This energy was presented by Buist et al. 28 and is given by e(q) = U g nHU + L g nHL + 1 2 HereH U =H U (A U (q 1 , U )) andH L =H L (A L (q 2 , L )) are general geometric terms representing the centers of mass of the upper and lower fluids, respectively (see Appendix A). From this energy, we can calculate the entropy variables as .
It was already proven in Buist et al. 28 that the flux terms of the TFM are energy-conserving, and the condition (15) takes entirely the same form for the PFTFM as for the TFM. Therefore, we state that (15) holds with This can be proven by simply evaluating the left-hand side (LHS) of (15), applying product and chain rules, and substituting geometric relations (A6) 2 . It is possible to substitute (A4) in (20), but we shall keep the current notation, in which the contributions from both fluids take the same form. Note that each term in (20) contains a multiplication with a velocity, so that h will be zero at closed boundaries, and there will indeed be no contribution to the global energy balance (18) in this case.

Determining the PFTFM free parameter based on energy conservation
The proof of energy conservation for the remaining terms in the PFTFM differs from the proof of energy conservation for the pressure terms of the TFM, since condition (16) differs from its TFM equivalent. Therefore, we reintroduce apostrophes to distinguish the PFTFM solution from the TFM solution, which may differ in caseK is set incorrectly.
This difference in solutions also stands in the way of proving energy conservation, but it turns out that these two issues can be resolved by a single approach. This approach is to setK based on the demand for global energy conservation. The LHS of (16) can be evaluated using ⟨v(q ′ ), d(q ′ )⟩ = Q ′ and substituting definitions (11), which yields ) .
For closed and inflow boundary conditions,Q ′ is determined by the boundary conditions in the same way for the PFTFM as for the TFM, so we haveK =Q and q ′ = q, and resubstituting (7) in (21) and by extension in the LHS of (16) leads to the conclusion that (16) holds with matching the TFM. However, due to arguments given in Section 2, for periodic boundary conditions generallyK ≠Q and q ′ ≠ q, and (7) cannot be substituted, since it is valid only for solutions of the TFM. In general, the RHS of (21) can be integrated over the domain to yield the contribution to the global energy balance (18): in which the advective and gravitational flux terms are encapsulated in the energy flux h (this was discussed above) and are therefore energy-conserving, and the expected contribution of the source terms will be examined in Section 3.4. Our key insight now is that (for periodic boundary conditions)K can be chosen to make the pressure terms of the model energy-conserving. To determine the value forK, we demand that the pressure terms are globally energy-conserving, meaning that they should not contribute to the global energy balance: Applying (4) leads toK .
This is the valueK should be set to in order to achieve global energy conservation, for periodic boundary conditions, wheṅ Q is unknown. We note that a similar expression was derived by Milewski et al., 25 for a two-layer shallow water model with a rigid lid, but the approximationK(t) = 0 was made for the analysis and computations. We will use the precise expression (25), meaning thatK is set as a function of the global solution at time t. With this choice, (23) reduces to (18). Using the proof that with this choice the PFTFM is equivalent to the TFM (see subsection 3.

3), (24) can be written as [Qp]
For the case of periodic boundaries, the new and improved PFTFM consists of (12) supplemented by (25), making it a combined system of differential and integral equations. For other boundary conditions, the PFTFM consists of (12) supplemented by a functionK =K(t), which is determined by the (possibly time-dependent) boundary conditions. These methods of settingK lead to a volumetric flow rate matching that of the original TFM solution:K(t) =Q(t) (see Section 3.3). It can be proven that with the changes made, and as long as consistent initial conditions are used, TFM solutions satisfy the PFTFM model equations, and vice versa (see Section 3.3). Therefore, the solutions to the PFTFM must be locally energy-conserving (with j = Qp) and momentum-conserving (see Appendix B).

Equivalence of PFTFM and TFM solutions
Given a solution q ′ (s, t), (25) provides a unique value forK(t) =Q ′ (t) in order to obtain a progression of the solution in time that is globally energy-conserving. In the derivation of the PFTFM from the TFM given in Section 2.2, the solution is unchanged, up to the step of setting the rate of change of volumetric flow rate as an explicit parameterK(t), instead of leaving it as a part of the solution in the formQ( q ′ ∕ t) (given by (8)). This will change the solution ifK ≠Q, but iḟ K =Q, we change nothing in this step, and the transformation from TFM to PFTFM can be freely reversed. Since we have found the unique expression forK =Q ′ that will progress the solution in an energy-conserving manner, and we know that the TFM is energy-conserving andQ will progress the solution q in an energy-conserving manner, it follows that, with identical initial conditions q ′ (s, t 0 ) = q(s, t 0 ), we haveK(t 0 ) =Q(t 0 ), and the progression of the solutions must be identical so that at all further times we have q ′ (s, t) = q(s, t) andK(t) =Q(t).
The equivalence of solutions to the TFM and PFTFM can be verified by checking that solutions to the TFM satisfy the PFTFM model equations, and vice versa. First consider a solution q(s, t), p(s, t),Q( q∕ t) =Q(t), which satisfies the TFM model equations (1) and the volume constraint (3), and its derived constraints. To check if this solution also satisfies the PFTFM model equations, we substitute q in the PFTFM system (12), and check if the following expression evaluates to zero: ) .
Clearly, if the following substitution can be made: then (26) reduces to the original TFM equations (which q and p were defined to satisfy), and therefore it evaluates to zero. For the TFM we know from (7) that (27) holds ifK is replaced byQ. In this case, it is trivial to setK(t) =Q(t), since we knowQ(t) as a part of the TFM solution. Therefore, the substitution given by (27) can be made, (26) reduces to zero, and the TFM solution satisfies the PFTFM model equations.
which satisfies the PFTFM model equations, and check if it satisfies the TFM model equations and its constraints. Substituting this solution in the TFM system (1) requires an expression for the pressure, so as a trial solution, we calculate the pressure based on (7) and the PFTFM solution. The result is: where we have used (13). This demonstrates that, using (7) to define the pressure gradient as a function of the solution, the PFTFM solution satisfies the TFM model equations, for anyK. However, this does not take into account the boundary conditions. For closed or mass inflow boundaries,K should agree with the inflow set at the boundaries. For periodic boundary conditions, p ′ should be periodic, meaning that leading again to the demand thatK should be set by (25). The proof is completed by referring back to Section 2.4, where it has already been demonstrated that the PFTFM satisfies the volume and volumetric flow constraints, as long as the initial conditions are consistent.
To summarize, ifK is set by (25), and the initial conditions satisfy the constraints, then any solution to the TFM satisfies the PFTFM, and any solution to the PFTFM satisfies the TFM. Therefore, we can conclude that the two models yield equivalent solutions.

Conservative formulation of the streamwise gravity source terms
A second addition to the PFTFM that we propose in this paper is an energy-conserving discretization of the source terms related to gravity in the streamwise direction. We first need to show that the continuous model conserves streamwise potential energy, before constructing a discretization to mimic this property (the latter will be discussed in Section 5.3). This addition extends upon the discretization presented by Buist et al. 28 and can be applied just as well to the TFM. The goal of this addition is to make the energy-conserving discretization more complete and extend its applicability to more general cases (that is, inclined flow). The source terms can be split into different components, one of which is the streamwise gravity term c g : The term c g can be incorporated in the energy-conserving formulation by adding a streamwise potential energy term to the energy (based on physical considerations): sin( (s ′ )) ds ′ , and dy ds = sin( (s)).

For this energy, the vector of entropy variables becomes
If e new is to be conserved, in the absence of other source terms, a new expression for h is required such that the following holds (compare to (15)): Regarding (16), nothing is changed, since only the terms in v that pertain to q 1 and q 2 are different than before, and (16) contains only terms pertaining to q 3 and q 4 . If (29) is satisfied, then the energy conservation equation becomes in which c rem is the remainder of the source terms; it contains all source terms except for c g . This equation differs from (17) through the changed expressions for e and h and a part of the source terms being removed from the right-hand side (RHS). It can be integrated over the domain to yield The extra terms on the LHS of (29), relative to (15), are given by Therefore, (29) is satisfied with and (30)  The derived h new is zero at closed boundaries, so that the boundary terms in (31) disappear and the global energy remains constant. However, h new is not periodic in the case of periodic boundaries with different y(s) at each boundary. The difference in y(s) between the boundaries makes the domain non-periodic. The change in energy resulting from the difference between h new (s 1 ) and h new (s 2 ) can potentially be balanced by a driving pressure gradient, which appears as a source term in the global energy conservation equation (31).
In conclusion, we have shown that the streamwise potential energy can be added to the energy, so that the amended energy is conserved (with appropriate boundary conditions). In contrast, remaining source terms such as the driving pressure gradient act as source terms in the energy equation.

Semi-discrete form of the TFM
We have shown how the continuous PFTFM, derived from the continuous TFM, can be adapted with an appropriate expression for the volumetric flow rate. This leads to solutions matching the TFM, and energy conservation. We will now use these insights to obtain a semi-discrete PFTFM with the same properties, making use of an existing energy-conserving discretization of the TFM. 28 To this end we must first derive the semi-discrete PFTFM from the semi-discrete TFM, mimicking the derivation of the continuous PFTFM from the continuous TFM in Section 2. We choose this route rather than directly discretizing the continuous PFTFM, because the chosen route ensures that the energy-conserving discretization of the TFM is applied to the PFTFM consistently (that is, that the correct discretizations of B and k are obtained). Such a consistent application of the TFM energy-conserving discretization to the PFTFM will ensure energy conservation, since this makes the semi-discrete PFTFM a reformulation of the energy-conserving semi-discrete TFM, for which energy conservation has already been demonstrated. 28 Therefore we first present the semi-discrete form of the TFM, as given in Buist et.al. 28 We have a staggered grid with N p 'pressure volumes' and N u 'velocity volumes', shown in Figure 2. Our semi-discrete vector of unknowns is based on its continuous equivalent (2), but includes an extra factor Δs so that the unknowns become quantities of mass or momentum. The quantities of mass are defined at the pressure volumes and the quantities of momentum are defined at the velocity volumes. With a line over the variable indicating interpolation: and with brackets indicating a jump: with The volume constraint closes the equations. Through differentiation to t and substitution of the mass conservation equations (with fluxes given by (35)), it implies the volumetric flow constraint with Just as in the continuous case, this constraint can be used to formulate an equation for the pressure. The discrete momentum equations can be summed as follows: Substituting (40) yields an equation for the pressure that compares to (7): Taking the difference of (42) with a shifted version of itself and applying (39) yields 1 Δs This is the discrete version of (9).

Semi-discrete form of the PFTFM
Having presented an energy-conserving discretization of the TFM and its constraints, we can derive the pressure-free version of it, following the approach of Section 2.2. Like in Section 2.2, we multiply (42) by d i to obtain with Substitution of (44) in (34) then yields the semi-discrete PFTFM where the apostrophes denote (dependence on) the PFTFM solution as opposed to the TFM solution (e.g. . HereK must be set as an input to the model, as discussed in Section 2.2.  Figure 2). At the boundary points, we prescribe the boundary conditions in weak form. If the first interior pressure grid point is given an index of 1 and the last is given an index of N p , this means that we prescribe dq 3,1∕2 ∕dt = 0, dq 4,1∕2 ∕dt = 0, dq 3,N p +1∕2 ∕dt = 0, dq 4,N p +1∕2 ∕dt = 0 (the initial conditions should be prescribed such that q 3,1∕2 = 0, q 4,1∕2 = 0, q 3,N p +1∕2 = 0, q 4,N p +1∕2 = 0 for closed boundaries). The boundary values of dq 1,1∕2 ∕dt, dq 2,1∕2 ∕dt, dq 1,N p +1∕2 ∕dt, dq 2,N p +1∕2 ∕dt then follow via an analysis of the characteristics, as explained by Sanderse and Veldman. 35 In case of inflow boundaries, the treatment is generally the same as for closed boundaries, but the time derivatives of the momenta can be set to some predefined function of time instead of zero (and similar for their initial values). The disadvantage of prescribing the boundary conditions in weak form as opposed to strong form is the introduction of a time integration error at the boundaries. However, weak boundary conditions are necessary for consistency with the prescription of the volumetric flow rate in weak form (through the parameterK). 27

Discrete aspects of enforcing the constraints
Regarding the connection between the volumetric flow rate and the boundary conditions, the same arguments apply as in the continuous case. Through definition (40), the volumetric flow rate is determined by closed or inflow boundary conditions, but not by periodic boundary conditions. Through (39) this is extended to the entire domain. For the TFM, the volume constraint (38) and the volumetric flow constraint (39) are enforced using a predictor-corrector algorithm involving the solution of the pressure Poisson equation (43), as a part of a semi-explicit Runge-Kutta method. 35 When running the TFM (for the purpose of comparing results to the PFTFM), we will solve the (implicit) pressure Poisson equation using a preconditioned conjugate gradient method. Satisfaction of the volume constraint is in turn necessary for total momentum conservation, and satisfaction of the volumetric flow constraint is necessary for energy conservation. 28 Total mass conservation is a property of the finite volume scheme and does not depend on the pressure Poisson equation.
For the PFTFM, just as in the continuous case, one can take the dot product of (46) with (1∕Δs)l and substitute (40) to findQ ′ i−1∕2 =K, and this holds for any index, which shows that the volumetric flow constraint (39) is satisfied naturally by the PFTFM, as long as the initial condition satisfies the constraint. Similarly, take the dot product of (46) with (1∕Δs) [ implying that the volume constraint is also naturally satisfied by the PFTFM, as long as the initial conditions satisfy the constraints. This shows how both models incorporate the volumetric flow constraint to determine the volumetric flow rate throughout the domain. The downside of the TFM is that it requires the solution of an extra equation (the pressure Poisson equation (43)), which adds to the computation time. Additionally, it introduces a numerical error which reduces the accuracy to which the constraints are satisfied, and as a result reduces the accuracy to which momentum and energy are conserved. In contrast, the time integration of the PFTFM can be performed in a fully explicit manner, and we will use a standard four-stage, fourth-order Runge-Kutta method. 38 The downside of the PFTFM is that the volumetric flow rate must be set explicitly throughK. This issue will be resolved for the discrete model in the following section, analogously to the continuous model.

Energy conservation for the semi-discrete model
To find an appropriate semi-discrete expression for the volumetric flow rate, we will use the demand for global energy conservation, as we have for the continuous model. Therefore we must study energy conservation for the semi-discrete model. Just as for the continuous model, a local energy conservation equation can be derived from the equations for mass and momentum conservation that constitute the model. The derivation starts by defining the discrete mechanical energy e i−1∕2 , from which the discrete vectors of entropy variables can be calculated: Taking the sum of the dot products of v i−1∕2,i−1 and v i−1∕2,i with the systems (46) for q i−1 and q i , respectively, yields the local energy conservation equation (with h i and j i as energy fluxes) under the following conditions: This means that if the LHS of (48) and (49) can be written in conservative form, that is, as the difference between an ingoing and an outgoing flux for the grid point i − 1∕2, then the scheme is energy-conserving. Therefore, (48) and (49) are the conditions for semi-discrete local energy conservation for the PFTFM, and are the discrete analogues of (15) and (16).
In case of periodic or closed boundaries, "integrating" (summing) (47) over the domain yields the global energy conservation equation We now consider a specific energy, and verify if the conditions for energy conservation are satisfied. Since the semi-discrete PFTFM is directly derived from the semi-discrete TFM, it should conserve the same energy (if the volumetric flow rate is set correctly). The energy which was shown to be conserved by the TFM in Buist et al. 28 is given by From this definition, the entropy variables can be calculated as Condition (48) takes exactly the same form for the PFTFM as for the TFM, and does not involve the pressure or the volumetric flow rate. Since our semi-discrete PFTFM was derived from the energy-conserving semi-discrete TFM that was shown to satisfy this condition, 28 and consequently uses the same expressions for the numerical flux (given by (35)), it is clear that this condition will be satisfied for the PFTFM as well. This can be verified via a direct evaluation of the LHS of (48), making use of a discrete product rule and a discrete counterpart of (A6). 3 This will show that (48) holds with Since the last terms in each line vanish for Δs → 0, this discrete definition of h is consistent with the continuous definition given by (20). Condition (48) being satisfied proves that the flux terms of the semi-discrete PFTFM are energyconserving.

Determining the PFTFM free parameter based on energy conservation
We now turn to condition (49), which needs to be satisfied for the "pressure terms" of the semi-discrete PFTFM to be energy-conserving. In contrast to condition (48), condition (49) is unique to the PFTFM, and it hinges on a correct expression forK. We again use apostrophes to distinguish the PFTFM solution from the TFM solution, and to distinguish functions of the PFTFM solution from functions of the TFM solution. With d i defined by (36), we obtain Therefore, similar to (21), the LHS of (49) evaluates to For boundary conditions for which the volumetric flow is known a priori (closed or inflow boundaries), andK is set accordingly so thatK =Q, we can resubstitute the discrete pressure equation (42) to show that (49) holds with j i = Q ′ p i , where we have used (39): the volumetric flow rate must be constant with s (and it is for the PFTFM, see subsection 4.4). Just as for the continuous model, in case of periodic boundary conditions,Q is an a priori unknown function of the solution. We again setK such that the considered terms conserve energy globally: Using (39) again, this leads toK At each time stepK should be set to the current value of the RHS of this expression. A standard explicit Runge-Kutta time integration method can be used, in which caseK should be calculated similarly at each stage. Setting the time derivative of the volumetric flow rate as described here corresponds to the 'option A' time discretization from Sanderse et al. 27 The alternative is to set K at each time step ('option B'), requiring a specific and careful adaptation of the Runge-Kutta method. This second option does not fit in our framework of semi-discrete energy conservation, of which the purpose is to ensure that the spatial discretization does not introduce an energy error. Therefore it is concerned only with how the different terms in the model relate to each other instantaneously at any given moment in time. In the semi-discrete model formulation onlyK is present, as opposed to K, and therefore we can only impose conditions onK.
The novel energy-conserving semi-discrete PFTFM consists of (46), with fluxes given by (35), and withK determined by (50) in case of periodic boundary conditions. For other boundary conditions,K is determined by the mass inflow set at the boundaries. Standard time integration methods can be applied to the semi-discrete PFTFM; we use a classic four-stage fourth-order Runge-Kutta method. 38 All the relations needed to implement the model are summarized in Appendix D.
With consistent initial conditions, the PFTFM will naturally satisfy the constraints (see Section 4.4) and conserve momentum exactly with time (see Appendix C). Via the same arguments as given for the continuous model in Section 3, the solution to the semi-discrete PFTFM will match that of the semi-discrete TFM, withK =Q. Between the fully discrete models, a very small difference remains, due to the details of the constraint-consistent time integration method of the TFM, involving a predictor-corrector algorithm and the solution of the pressure Poisson equation (43), instead of a direct evaluation of (50).

Conservative discretization of the streamwise gravity source terms
In Section 3.4 it was shown that for the continuous model, the streamwise potential energy can be included in the local energy conservation equation. Only a carefully chosen discretization will replicate this property in the discrete setting. In this subsection, we shall pose such a discretization and show that it conserves the streamwise potential energy together with the normal potential energy and the kinetic energy. We propose the following discretization for the streamwise gravity terms: implying that the channel inclination has been discretized as The streamwise potential energy is added to the local energy, which becomes Other discretizations for the source terms and the energy are conceivable, but we will show that this particular set is energy-conserving, in conjunction with the discretization of the core model that was shown to be energy-conserving in the previous subsection.
Our chosen discretization implies that y is defined on the pressure grid, and sin( ) is defined on the velocity grid. This choice is natural since it implies that the streamwise potential energy is defined on the pressure grid and is interpolated as a whole to the velocity grid, which matches the treatment of the normal potential energy terms. This is particularly useful at the boundaries in a nonperiodic domain, where velocity half-volumes are present which require this interpolation to proceed differently (see Figure 2). The potential energy of these half-volumes is determined only from the nearest interior grid point: e p,1∕2 = (1∕2)e p,1 and e p,N p +1∕2 = (1∕2)e p,N p . Treating the boundaries this way ensures that the act of interpolating the potential energy from the pressure to the velocity grid does not change the global potential energy. The kinetic energy is calculated directly from the boundary values without interpolation (masses and momenta are both available here).
In the proof of energy conservation, condition (48) is amended to The discrete condition on j, given by (49), remains unchanged. If the conditions are satisfied, the energy conservation equation becomes in which c rem,i contains the remaining source terms. This equation can be integrated over the domain to yield the new global energy conservation equation in which the streamwise gravity terms have been removed from the right-hand side. For the energy given by (51), the entropy variables are Evaluating the LHS of (52) shows that there are only a few additional terms compared to those that were already shown to be conservative in the previous subsection. It can be shown that these additional terms can be written as Therefore (52) is satisfied with which remains consistent with the continuous expression (32), and completes the amended local energy conservation equation (53).
In conclusion, we have shown that with a particular discretization of the source terms related to streamwise gravity, and a particular definition of the discrete streamwise potential energy, it can be shown that this energy is conserved in conjunction with the normal potential energy and the kinetic energy. Therefore, with this particular discretization, the conservation properties of the continuous model (which were discussed in Section 3) apply to the semi-discrete model.

NUMERICAL EXPERIMENTS
We conduct numerical experiments in order to confirm the properties of the proposed discretization of the PFTFM. The first test case has periodic boundaries. It will be used to test the new expression forK and its implementation in a fully explicit Runge-Kutta method. The test will demonstrate that the TFM and PFTFM solutions are now equivalent up to a small numerical error. Additionally, it will show that with the proposed adaptation, momentum, and energy are conserved up to high precision. The second test case has closed boundaries and an inclined domain, and will be used to test the proposed conservative discretization of the streamwise gravity terms. Both test cases are chosen such that for the continuous model the global energy remains constant, so that the discretization can be tested by verifying that The spatial discretization should conserve energy exactly, with a residual temporal error that decreases as the time step is reduced. The nondimensional energy error will be defined as h the energy of the initial condition. The test cases are limited to the 2D channel geometry, since the discretization is not exactly energy-conserving for the circular pipe geometry (see the footnote in Section 5.1). The parameters of the cases are chosen such that they avoid discontinuities, which would require the addition of diffusion to the numerical scheme in order to mitigate oscillations.
The TFM with the hydrostatic pressure assumption, as given in Section 2.1, is only conditionally hyperbolic. 10,33 This means that, for certain configurations of the flow variables, the eigenvalues can become complex. As a result, a linear stability analysis will show an unbounded growth rate for perturbations of vanishing wavelength, which is unphysical, and can lead to a lack of convergence of numerical solutions. 11,39 Since the PFTFM has the same eigenvalues as the TFM, 27 it retains this property. In the test cases considered here, the region of state space where the model attains complex eigenvalues is avoided by the choice of parameters and initial conditions. This implies that the velocity difference between the lower and upper fluid is kept small (below the inviscid Kelvin-Helmholtz stability criterion).

Traveling wave
This case is a modified version of a case featured in Buist et al. 28 In this case, we initialize a sine-shaped perturbation to two fluids in steady state stratified flow in a (level) periodic domain. The parameters of the case are inspired by the Thorpe experiment, 40 and are given in Table 1.
these perturbations give rise to exact solutions to the linearized model 11 of the form where the relative magnitudes of the perturbations are given by

Δ̂L.
We set Δ̂L = 0.05 and calculate the remaining components of the perturbation according to the expressions given above, selecting the value for ∕k corresponding to the plus sign in (54). A projection step is performed on this perturbed state, to yield a set of initial conditions that satisfies the constraints. In the linear approximation, these initial conditions lead to a single wave traveling to the right at constant amplitude. In the full nonlinear model, the wave deforms over time, which can be seen by plotting the solution at integer multiples of the wave period Lk∕ (Figure 3). Figure 4 shows the difference between the TFM and the PFTFM results for two components of the global solution: the potential energy and the volumetric flow rate. The results show that by settingK as described in Section 5.2, the correct volumetric flow rate is obtained by the PFTFM, which matches that of the TFM (for reference, the initial value of the volumetric flow rate is Q 0 = 3.3 ⋅ 10 −2 m 2 s −1 , and the variation in Q over the course of the 30 second simulation is of order

F I G U R E 5
Conserved quantities for the traveling wave test case. Left: potential, kinetic and total energy relative to their initial values.

Right:
( can be viewed at wileyonlinelibrary.com] 10 −7 m 2 s −1 ). The agreement in potential energy is also excellent. This resolves the issue reported in Sanderse et al., 27 where the authors observed an unjustified difference in volumetric flow rate between TFM and PFTFM solutions, with an accompanying difference in local solutions. Having repaired the volumetric flow rate, the simulations exhibit near machine precision agreement in local solutions between TFM and PFTFM (not shown). This agreement is independent of time step. This near-identical solution is obtained at a significant reduction of computational expense: the wall-clock time for the PFTFM is 232 s while for the TFM it is 733 s. This is roughly consistent with the results shown in Sanderse et al., 27 where it was already concluded that the PFTFM saves computing time by not having to solve the Poisson equation for the pressure. This conclusion proves to remain valid even though in that work a direct solver was used for the pressure Poisson equation instead of the iterative (conjugate gradient) solver used here.
The energy behavior is shown in Figure 5, for a simulation with a resolution of N p = 40 and a time step of Δt = 5 ⋅ 10 −4 s. The exact solution of the linearized model has constant kinetic and potential energy since it is an unchanging traveling wave, but the computed solution shows a somewhat erratic exchange between kinetic energy and a potential solution, at a period of a bit under four times the wave period. This exchange is small relative to the base energy, around the level of the spatial numerical error. Despite this irregular exchange, the total energy is conserved to great precision, confirming the energy-conserving properties of our discretization. IfK is not set as described in Section 5.2, but is set incorrectly toK = 0, this leads to a total energy error that is many orders of magnitude larger than the error shown here.
It was discussed in Section 4.4 that the absence of the Poisson equation for the pressure in the formulation of the discrete PFTFM may lead to improved momentum and energy conservation when compared to the TFM. This is demonstrated in Figure 6, where we show the conservation errors made by both models, varying the time step. Both models are solved with a resolution of N p = 40, and for the Poisson equation the TFM uses a preconditioned conjugate gradient solver set to a tolerance at the order of machine precision (10 −15 ), so that it converges to the smallest possible residual (within a maximum of 50 iteration steps). The error in momentum conservation is constant with time step for both models, but is consistently larger for the TFM, and this can be ascribed to an increased volume constraint error. Regarding the energy, both models initially converge to fourth order with time step, but the PFTFM converges further, before also levelling off at a significantly smaller error than the TFM. This indicates that the remaining error in the pressure Poisson solve forms a limiting factor to the accuracy of energy conservation for the TFM.

Sloshing in a closed tank
This is a modified version of a test case presented by Sanderse and Veldman. 35 In this test case, a section of duct is bounded at the left and right ends by closed boundaries, forming a closed tank. Two fluids, with the density of water and air, are initialized at rest with a horizontal boundary between them. At t = 0, the tank is set at an incline, causing the water to start flowing towards the lower end of the tank, pushing the air toward the higher end. This case tests the performance of our discretization of the source terms associated with the streamwise gravity, which was given in Section 5.3. Here it was shown that with an appropriate discretization of these terms, and including the streamwise potential energy in the energy definition, the streamwise gravity does not act as a source term in the energy balance. Therefore, we expect to find spatially exact conservation of energy.
The parameters of this case are given in Table 2. The experiment differs from the reference case 35 due to the absence of wall and interface friction, the reduction of the inclination from 2 to 1 degrees, and the change of duct geometry from circular pipe to 2D channel. The number of pressure cells taken for the discretization is N p = 80, just as in the reference case. 35 Figure 7 shows the different components of the solution in a space-time plot. At the start of the simulation, the results show waves in the hold-up emanating from both sides of the tank. From the right boundary, an expansion wave travels to the left, reducing the hold-up. From the left boundary, a compression wave travels to the right. The expansion wave is faster than the compression wave, reaching the left boundary around t = 1.3, while the compression wave reaches the right boundary at t = 1.7. Figure 8 shows the solution at the points in time when a wave front meets a boundary. The waves reflect at the boundaries, and when they meet again the compression wave gains amplitude. This repeats itself at future encounters between the two waves, building up toward a shock wave. The discontinuity which is being formed leads to the onset of numerical oscillations, which can be observed starting around t = 4 s. The numerical oscillations may be eliminated through the addition of physical (molecular and turbulent) viscosity, 39 (grid-independent) artificial diffusion, 19,42 or (grid-dependent) numerical viscosity. 37 This is outside the scope of this paper, since here we are interested in demonstrating the energy conservation property of the discretization, which would be spoiled by adding viscosity. Our discretization can serve as an energy-conserving basis to which energy-dissipating terms may be added. Figure 9 shows that, with a small time step of Δt = 10 −4 s, the total energy is conserved up to machine precision while potential and kinetic energy are exchanged. The potential energy includes the normal potential energy and the streamwise potential energy. It is due to the streamwise potential energy present in the initial condition, which is converted into kinetic energy, that motion is initiated. In these simulations there is no damping (except for minute damping due to time integration error), so the fluids will continue sloshing for as long as the simulation is run. Figure 10 shows that the energy error converges to high precision with time step, while remaining roughly constant with grid resolution. The convergence with respect to time is approximately fourth order, corresponding to the order of the Runge-Kutta time integration scheme. This demonstrates that the spatial discretization is exactly energy conserving and that the remaining error is only temporal.

F I G U R E 9
Conserved quantities for the sloshing test case. Left: potential, kinetic and total energy relative to their initial values. Right: Finally, we verify that the solution of the PFTFM matches the solution of the TFM. Since the volumetric flow rate is known a priori (Q = 0 andQ = 0), this is achieved without difficulty. This is demonstrated in Figure 11 by comparing the potential and kinetic energy obtained by the TFM and the PFTFM. The difference between the two models is near the order of machine precision (and the same holds for the local solutions). Again the PFTFM saves computing time by not having to solve the Poisson equation for the pressure: the wall-clock time for the PFTFM is 208 s while for the TFM it is 665 s.

CONCLUSIONS
In this article, we have removed the discrepancy between the PFTFM and the original two-fluid model from which it is derived. For periodic boundaries, this requires setting the rate of change of the volumetric flow rate based on the demand of global energy and momentum conservation. This change resolves the issues of disagreement between PFTFM and TFM solutions and of energy conservation in a combined approach. Our approach may inspire similar improvements to other pressure-eliminated versions of the TFM, such as the FFM. The PFTFM stands out from other pressure-eliminated versions of the TFM in that mass, momentum, and energy-conserving discretizations of the TFM can be directly applied, and in the fact that it preserves the correct shock relations. A further addition to the energy-conserving framework of Buist et al. 28 was made: the streamwise potential energy and its related source terms were included in the energy-conserving formulation. The amended energy is conserved in inclined domains without friction. We have presented an energy-conserving discretization of these source terms and their associated energy, extending the energy conservation properties of the previously presented discretization.
It was shown how a semi-discrete version of the PFTFM can be derived from the TFM, directly demonstrating how the energy-conserving discretization of the TFM presented in Buist et al. 28 can be applied to the PFTFM. In the framework of this semi-discrete PFTFM, a discrete equivalent of the new expression for the (rate of change of the) volumetric flow rate was derived, based on the demand for discrete global energy conservation. This addition completes the energy-conserving discretization of the PFTFM for the important case of periodic boundary conditions. With the proposed discretization, the discrepancy between PFTFM and TFM solutions, observed in numerical simulations in Sanderse et al., 27 is removed for the discrete model.
The PFTFM modified and discretized as shown in this work conserves mass, momentum, and energy, without needing to account for implicit constraints, making the time integration fully explicit (standard explicit Runge-Kutta methods can be applied). Thus solutions are obtained that are equivalent to solutions of the TFM (down to machine precision), at significantly reduced computational cost. Removing the error attributable to the pressure Poisson solve in the TFM leads to slightly more precise conservation of momentum and energy for the PFTFM. A small temporal error remains present in the energy, which may be resolved in future work by developing an energy-conserving time integration method for the PFTFM.
In their current form, the model and its discretization contain no numerical or physical diffusion. This omission is necessary to achieve exact conservation of energy, but has the effect that discontinuities cannot be correctly resolved. To resolve this, our model could be extended with diffusion terms as found in, for example, Fullmer et al., 39 by simply including them in the source term vector. These should be a strictly dissipative addition to the energy-conserving discrete model, and should produce the dissipation required at discontinuities, without excessive dissipation where the solution is smooth.
The elimination of constraints and of the nonconservative pressure terms brings the model closer in mathematical form to related models such as the TLSWE. Like the TLSWE, the TFM is only conditionally hyperbolic. This can be explained by the fact that both models describe the flow of two discrete fluids with a sharp interface, under hydrostatic balance. The PFTFM and the modifications we propose in this work do not resolve the hyperbolicity issue. However, the PFTFM's energy-conserving formulation, and its increased mathematical similarity to the TLSWE, may facilitate the application of techniques developed for dealing with complex eigenvalues and instability for the TLSWE. 43,44

ACKNOWLEDGMENTS
This work was supported by the research program Shell-NWO/FOM Computational Sciences for Energy Research (CSER), project number 15CSER17, which is partly financed by the Netherlands Organization for Scientific Research (NWO).

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. 1 To see this, note that this model can also be obtained by taking the product of A and (1). A is a projection matrix satisfying A 2 = A, and it is singular. The term A q∕ t needs to be split in the following manner: A q∕ t = q∕ t − B q∕ t = q ′ ∕ t − kK. 2 In this derivation g n is assumed to be constant along s. 3 Unlike their continuous versions, the discrete geometric conditions are only satisfied exactly by geometries with d 2 H L ∕dA 2 L = 0 (for example the 2D channel geometry). For other geometries, such as the cylindrical pipe geometry, these conditions are satisfied up to a numerical error of order (Δs) 3 . 28 4 Using the Churchill friction model 41 with viscosities of U = 1.5 ⋅ 10 −3 kg m −1 s −1 and L = 1 ⋅ 10 −3 kg m −1 s −1 . 5 For the channel geometry, we substitute A = H and P int = 1.

APPENDIX A. GEOMETRIC RELATIONS
The equations are written in terms of the cross-sectional areas occupied by each fluid, which in general can be defined to be related to the interface height H L via with w(h) the local duct width. Note that w(H L ) = P int , where P int is the (generalized) interface perimeter which is shown for a circular pipe geometry in Figure 1 where we have substituted A L = A − A U . For the pipe geometry, we make the transformation h = R (1 − cos( * )), with * the integration variable and the wetted angle, to get 34 BesidesĤ L andĤ U , the following geometric quantities are used in the energy definition: We can apply Leibniz' rule to (A1) to obtain where we have used H L + H U = H. Applying the same technique to (A2) yields A similar calculation, in which we assume that H is constant, yieldŝ Finally, comparison to (A4) yieldsĤ

APPENDIX C. MOMENTUM CONSERVATION FOR THE DISCRETE PFTFM
Just as in the continuous case, it can be shown that our expression (50) forK (applied in the case of periodic boundary conditions) leads to total momentum conservation (in addition to energy conservation). The semi-discrete equation for the conservation of local momentum m ′ i−1∕2 = m(q ′ i ) = q ′ 3,i−1∕2 (t) + q ′ 4,i−1∕2 (t) can be derived by taking the inner product of z = [ 0 0 1 1 ] and (46): and the global momentum equation is derived by summing over the domain: Here we have substituted ⟨ z, d ′ i ⟩ = A U,i−1∕2 + A L,i−1∕2 = A, which means that the volume constraint (38) must remain satisfied, as it is for the PFTFM (see Section 4.4). For periodic boundary conditions, settingK according to (50) leads to global momentum conservation, in absence of source terms. WithK =Q, substituting (42) leads to which is a momentum balance matching that of the TFM. In contrast to energy, total momentum is conserved regardless of the discretization of the numerical fluxes, due to the nature of the finite volume scheme (just like total mass). In addition, it is conserved by a simple forward Euler time discretization, and by extension also by explicit Runge-Kutta methods. The proof for this is simple since z contains only constants, so that the following holds:

APPENDIX D. SUMMARY OF THE DISCRETE MODEL
We summarize all the relations needed to implement the energy-conserving PFTFM. The basis of the model is (46), which describes the evolution of the variables (33), which are defined on the staggered grid pictured in Figure 2. The matrices A i and B i and the vector k i are defined by the vectors d i and l T according to (45), which are in turn defined by (36) and (41). The fluxes f i are defined by (35), and the difference of these fluxes should be taken: The source term vector, including a driving pressure gradient and streamwise gravity, is given by (37). The PFTFM does not include constraints or an equation for the pressure. Instead it is necessary to define the rate of change of the volumetric flow rate,K. For closed boundaries, this is simplyK = 0. For an inflow boundary, we havė where q 3,1∕2 and q 4,1∕2 are the (prescribed) left boundary values of q 3 = U u U A U and q 4 = L u L A L , respectively. For periodic boundaries,K is given by (50). For closed or inflow boundaries, the boundary values of dq 3 ∕dt and dq 4 ∕dt are prescribed, and the boundary values of dq 1 ∕dt and dq 2 ∕dt follow from a characteristic boundary treatment described in Sanderse and Veldman. 35 The system (46) describes the problem locally and should be used to compose the global problem, over the whole domain. The resulting system, consisting of the semi-discrete governing equations with the semi-discrete expression foṙ K, can be integrated in time using standard explicit time integration methods. In this work, we use a classic four-stage fourth-order Runge-Kutta method. 38 Further details on the time integration, including how to prevent constraint drift, are given in Sanderse et al. 27