Energy-conserving formulation of the two-fluid model for incompressible two-phase flow in channels and pipes

We show that the one-dimensional (1D) two-fluid model (TFM) for stratified flow in channels and pipes (in its incompressible, isothermal form) satisfies an energy conservation equation, which arises naturally from the mass and momentum conservation equations that constitute the model. This result extends upon earlier work on the shallow water equations (SWE), with the important difference that we include non-conservative pressure terms in the analysis, and that we propose a formulation that holds for ducts with an arbitrary cross-sectional shape, with the 2D channel and circular pipe geometries as special cases. The second novel result of this work is the formulation of a finite volume scheme for the TFM that satisfies a discrete form of the continuous energy equation. This discretization is derived in a manner that runs parallel to the continuous analysis. Due to the non-conservative pressure terms it is essential to employ a staggered grid, which requires careful consideration in defining the discrete energy and energy fluxes, and the relations between them and the discrete model. Numerical simulations confirm that the discrete energy is conserved.


Introduction
The one-dimensional (1D) two-uid model (TFM) is a dynamic model for strati ed ow in channels and pipes.It simpli es the full three-dimensional multiphase ow problem by resolving only the cross-sectionally averaged quantities (hold-ups, velocities, and pressure), which are often of practical interest.There are many variants of the model, but the basic idea, of two interacting uids whose behaviour is cross-sectionally averaged to obtain a 1D model, was introduced by Wallis (1969) [34] and Ishii (1975) [14].The model has among others applications in the oil and gas industry [12], in CO 2 transport and storage [3], and in nuclear reactor safety analysis [4].
An unsolved issue with the basic version of the TFM is that the initial value problem for the governing equations is only conditionally well-posed [22].This means that it is well-posed for some ow con gurations and ill-posed for others (e.g. when there is a large velocity di erence between the two uids).Conventionally, ill-posedness of the TFM is demonstrated by a linear stability analysis which shows an unbounded growth rate for the smallest wavelengths, when the values of the model variables are such that the eigenvalues are complex.In this case the solution is said to carry no physical meaning [19].However, when drawing conclusions on the well-posedness of the TFM, it is important to also consider its nonlinear aspects, and not only rely on a linearized analysis [17,29].Examples of studies that have included nonlinear e ects in the TFM analysis can be found in [16,20].However, a complete nonlinear analysis, with implications for obtaining a robust discretization, is still missing.
In this work, we strive towards such a nonlinear analysis by presenting an expression for an energy which is conserved by the full (nonlinear) TFM, in its incompressible and isothermal form.This approach is motivated by the fact that for the incompressible Navier-Stokes equations such an analysis provides stability estimates [8,26], and that for compressible equations it is closely related to the concept of entropy stability [30].Important to note is that such an energy is not the thermodynamic energy for which a separate conservation equation exists in the compressible TFM.Rather, the considered energy conservation is an inherent property of the mass and momentum conservation equations that constitute the incompressible TFM: the energy is a secondary conserved quantity of the model.Its physical meaning is therefore the mechanical energy of the system (kinetic plus potential energy).
In order to derive this mechanical energy equation, we take the approach from [11], in which the dot product of the shallow water equations (SWE) and a vector of entropy variables is taken in such a way that a scalar energy equation results.However, an important di erence with the SWE (and two-layer SWE [10]) is the presence of non-conservative pressure terms that are linked to the constraint that the uid phases have to ll the cross section.Another important di erence is that we consider arbitrary duct geometries, as opposed to the 1D SWE which in e ect utilizes a planar channel geometry.Given these di erences, the key challenge is thus to nd a conserved energy and corresponding energy ux function for the TFM, and this will be the rst main focus of this paper.
The second focus of this paper is to derive a spatial discretization which conserves a discrete version of the energy.Again, our approach is inspired by methods which have been developed for the SWE [9].An important di erence is that these methods are designed for collocated grids, while we will adapt them to a staggered grid.This is motivated by the presence of the (non-conservative) pressure terms in the TFM, which makes the use of a staggered grid much more convenient (similar to the case of the incompressible Navier-Stokes equations [8]).However, the staggered grid introduces new challenges, for example in terms of the de nitions of the energy and energy uxes.We will derive a discretization method that tackles these issues and propose a new set of numerical uxes on a staggered grid that are energy conservative.This discretization can also be viewed as an extension of the SWE discretization found in [33], where a di erent method is used to obtain a mass-, momentum-, and energy-conserving discretization on a staggered grid.
This paper is set up as follows.First, in section 2 we present the governing equations of the TFM.In section 3 we discuss the conditions for energy conservation, and introduce an energy and energy ux that satis es these conditions, providing local and global energy conservation equations for the continuous TFM.We outline how the equations are discretized in section 4, while leaving open the speci c form of the numerical uxes.Then, in section 5, we present the discrete versions of the continuous conditions for energy conservation, and propose a set of new conservative numerical uxes.Finally, in section 6 we present numerical results which demonstrate exact conservation of the aforementioned energy.

Governing equations
The 1D TFM, as considered in this work, describes the separated ow of a (heavier) lower uid and a (lighter) upper uid through a channel or pipe.It can be derived by applying a cross-sectional averaging procedure to the Navier-Stokes equations [15,29].An important assumption made in the derivation of the model is that the streamwise length scale is much larger than the normal length scale (i.e. the pipe diameter), which is referred to as the long wavelength assumption.As a consequence, along the normal direction the ow is in hydrostatic balance.We will omit source terms, such as wall friction, since such terms are sources or sinks of energy, and we are interested in the energy conservation properties of the core model.Good discussions of the assumptions underlying the TFM are given by [23,24].The cross-sectionally averaged equations can be written in the following concise form [27,28] (with = ( , )): where constitutes the vector of 'conserved' variables 1 , namely the mass and momentum of each phase: Here and are the densities, and are the cross sections, and and are the averaged velocities, all of the upper and lower uids, respectively.We consider the isothermal, incompressible case, so that and are constant.
The uxes describe convection of mass and momentum and gradients in the interface level.In terms of they are given by where = ( ) and = ( ) are geometric terms (to be discussed shortly), and is the gravitational acceleration in the normal direction.
The fth variable is the interface pressure , and the non-conservative pressure terms are given by ( ∕ ) with ( ) = 0 0 The quantities = ( ( 1 , )) and = ( ( 2 , )) are geometry-dependent and are de ned by Here the di erence between the coordinate ℎ and the two-uid interface height is integrated over the area occupied by the upper uid and the area occupied by the lower uid, respectively.Using these general expressions, the model equations are valid for arbitrarily shaped cross sections.See Appendix A for evaluations of the integrals for the 2D channel and circular pipe geometries.The spatial derivatives of and that appear in the uxes are known as the level gradient terms, which result from the hydrostatic variation of the pressure.
Since the upper and lower uid together ll the pipe, the system is subject to the volume constraint The entire system therefore consists of four evolution equations plus one constraint, and four 'conserved' variables plus the pressure.In our incompressible setting, a derived constraint can be obtained by di erentiating the constraint (6) and substituting the mass equations, leading to [28]: which can be integrated in space to give that the volumetric ow is constant in space, and a function of time only: This derived constraint, termed the volumetric ow constraint, can be seen as the incompressibility constraint for the TFM.We can use these constraints to set up an equation for the pressure.The pressure equation is obtained by summing the momentum equations [28]: which can be expanded and rewritten with the de nition of to yield Finally, taking the derivative of this equation to and applying constraint (7) gives This is a 'Poisson-type' equation for the pressure, which can be used in place of (6) to close the system of equations.In our numerical algorithm (discussed in section 4) we apply a discrete version of (11) in this manner.

Outline: conditions for energy conservation
Having set-up the TFM governing equations, the rst key objective of this paper is to prove local and global energy equalities that are implied by this equation set.This is similar to the energy analyses for e.g. the incompressible Navier-Stokes equations [8], the SWE [31], and the two-layer SWE [10].In all these models, no energy conservation equation is included in the model, but energy conservation follows from the mass and momentum conservation equations alone.It can therefore be said that the energy is a secondary conserved quantity.
Our proof of global energy conservation follows the approach in [9,11] and starts by showing that a local energy conservation equation of the form can be derived, purely based on manipulating the governing equations (1).Here ( ) is the local energy, and ℎ( ) and ( ) are energy uxes (to be detailed later); ℎ( ) is not be confused with the normal coordinate ℎ shown in Figure 1.If ( 12) holds, then it can be integrated in space to yield where the last equality ('=0') holds in case of periodic or closed boundaries, and the global energy ( ) is de ned as The key is therefore to obtain the local energy conservation equation (12).To achieve this, one rst postulates an energy ( ) (typically guided by physical considerations).Second, one calculates the vector of so-called entropy variables, de ned as2 ( ) ∶= .
Taking the dot product of the system (1) with leads to in which we have ignored source terms (as indicated before), and the brackets denote a dot product over the vector elements: ⟨ , ⟩ ∶= .
The time derivative term can be written as so ( 15) becomes an equation for the time evolution of the energy.
Given an expression for , the art is to nd an energy ux ℎ that satis es since then the second term in (15) can be written in the (locally) conservative form given by (12).In order to get a condition solely referring to the relations between di erent functions of (i.e.independent of ), the chain rule (valid for smooth solutions) is employed to convert (17) to: This is the condition encountered in e.g.[9] and [11] for an energy ux ℎ to conserve a given energy (or, more generally: entropy function) of the SWE.Likewise, we need to nd a ux such that the product of and the pressure gradient can be written in conservative form: The di erence between ℎ and lies in the fact that ℎ is responsible for the spatially conservative terms of the governing equations, whereas takes the non-conservative part into account.Perhaps surprisingly, we will show that these non-conservative terms ( ∕ ) can indeed be written in conservative form in the energy equation.An alternative formulation of condition ( 19) is given by In order for the local energy to be conserved, there must exist a (for the given and resulting ) such that this condition is satis ed.An important di erence between this derivation and the derivation for the SWE as found in e.g.[11] is the non-conservative pressure term.Although the two-layer SWE [1] also features a non-conservative term, in the TFM the non-conservative term depends on a variable for which there is no evolution equation (namely the pressure).This pressure term is instead linked directly to the volume constraint (6) and volumetric ow constraint (8) [28], which are not present in the SWE.For a system in conservative form without source terms, (18) is the only condition.This condition is emphasized in literature (e.g.[18]) as the condition for the existence of an entropy function.The derivation of energy conservation for the conservative part of the TFM system thus matches the derivation of an entropy condition for a conservative hyperbolic system.
In summary, the task is to nd a set , ℎ and, which satisfy conditions (17) and (19) for the current model with ux and pressure terms ( ∕ ).The alternative conditions (18) and (20) yield results more directly and will therefore be used in the following section.The result is the local energy conservation equation (12), and global energy conservation then follows directly.

Choice of energy and energy uxes
We will show that the energy is conserved by the TFM (in absence of source terms).Here = ( ( 1 , )) represents the center of mass of the upper uid multiplied by and = ( ( 2 , )) represents the center of mass of the lower uid multiplied by (see Appendix A), so that the rst two terms can be recognized as the potential energy of the upper and lower uid, respectively.The third and fourth terms represent the kinetic energy of the upper and lower uid, respectively.Therefore, this energy has a clear physical interpretation.
The entropy variables are given by ( ( 1 , )) and = ( ( 2 , )) representing the uid layer thickness of the upper and lower uids, respectively (see Appendix A).It is important that in the energy and energy ux terms concerning the upper uid we use ( ), ( ), and ( ), while for the lower uid we use ( ), ( ), and ( ).It is possible to use the volume constraint to change this functional dependence, but our choice leads to an elegant form of the energy conservation conditions.
The task is to nd ℎ and .We start with : the pressure term needs to satisfy (20).Straightforward evaluation gives ⟨ , ⟩ = with the volumetric ow rate given by (8).Because of the volumetric ow constraint (8), the second term of (20) vanishes, so that the condition on the pressure gradient evaluates to and ( 20) is satis ed with = . ( We note that is the pressure that enforces incompressibility -it does not include a driving pressure gradient (which would appear as a source term in the TFM governing equations).Therefore is periodic in space in the case of periodic boundaries.In the case of closed boundaries, must be zero, meaning that = 0 throughout the domain.This means that when integrating (12) over a closed or periodic domain, the terms involving vanish, and thus this de nition for is compatible with global energy conservation as described by (13).
The next task is to nd ℎ.Based on the form of ℎ for the SWE and condition (18), we propose the following choice which can be shown to satisfy condition (18) by computing: The last two entries in these vectors are equal because of relations (A.10), derived in Appendix A. The rst two entries are equal due to the geometric relations (A.6), which we repeat here in terms of the conserved variables : These relations follow directly from the de nitions of these geometric quantities and hold for arbitrary duct geometries.Note that, alternatively, condition (17) can be used (instead of ( 18)), which leads to the following conditions: which may also be shown to be satis ed directly via application of Leibniz' rule to the de nitions of and .These last two conditions will play an important role in the discrete analysis in section 5.
In conclusion, we have proposed a novel set of , ℎ, and for the TFM and have shown that the local energy conservation equation ( 12) is satis ed.

Reformulation in terms of the entropy potential and conditions on uxes
Conditions ( 17) and (19), or their alternatives ( 18) and ( 20), were used in the previous section to nd a combination of , ℎ and for the continuous TFM, given the uxes from the governing equations.In section 5, we will instead aim to nd discrete ux functions, given discretized versions of , ℎ and (that are inspired by their continuous counterparts).Equation ( 17) is not a very useful formulation to nd such numerical ux functions, because it is a condition imposed on the jump in , rather than itself.Therefore, equation ( 17) is reformulated using the concept of the entropy potential [9,31].
The entropy potential is de ned to be related to , , and ℎ in the following manner: With this de nition, we can reformulate condition (17) using the product rule ( ⟨ , ⟩ = ⟨ , ⟩ + ⟨ , ⟩) The entropy potential can be directly calculated from its de nition (27) and is given by: Because this entropy potential is based on an ℎ that satis es ( 17), ( 28) is satis ed by construction.Nevertheless, we outline the details to convert (28) into conditions on the individual numerical uxes, since they will be exactly mimicked by our discrete analysis in section 5. We rst introduce the following notation for the uxes, and split them into the following components: Here 3, and 4, are the momentum advection terms, and 3, and 4, are the level gradient terms (divided by ).These uxes and the de nitions for ( 22) and ( 29) can be substituted in (28).The resulting condition is rst split into two conditions: one condition proportional to , and one not proportional to .This is done on the basis that the mass and momentum advection terms do not depend on in the continuous case (see ( 3)), and should not depend on in the discrete case.These two conditions are split again on the basis that 1 and 3 should not depend on 2 and 4 , and 2 and 4 should not depend on 1 and 3 .We obtain the following four conditions: As mentioned, these equations are by construction satis ed by the ux vector (3).One important remark is that after we reformulated in terms of , the geometric conditions (26) encountered in subsection 3.2 still need to be satis ed in order for (31c) and (31d) to hold.

Comparison of the energy and energy uxes to those of other models
Here we compare the expressions obtained for and ℎ to results from literature for other models, focusing on the case of a channel geometry.The expression (21) for for the channel geometry can be obtained by substitution of the channel-speci c evaluations of and (Appendix A): For a single layer uid, such as the single layer SWE, only the third and fth terms remain, and they are consistent with the SWE entropy function as discussed in [11] (without channel inclination).
To compare with two layer SWE theory, we rewrite (32) using the volume constraint ( 6) to obtain This is the expression presented by Abgrall and Karni (AK) [1] and Fjordholm [10] as an entropy function for the two-layer SWE.The energy found in the present study can be seen as a generalization of the two-layer SWE energy to arbitrary duct geometries.When comparing our energy ux ℎ for the TFM to the one for the two-layer SWE, it should be realized that the two-layer SWE can be obtained from the TFM by the choice = .This means that the pressure ux of the TFM needs to be added to ℎ in order to compare with the SWE expressions.In our notation, the two-layer SWE entropy ux given by [1] is Our expression for ℎ for a channel is given by Upon adding = = to ℎ ch , and after some rewriting, we see that our TFM energy ux is consistent with the two-layer SWE entropy ux: To conclude, our proposed energy (21) and energy ux ( 24) can be seen as a generalization of the two-layer SWE energy and energy ux to arbitrary duct cross sections.

Semi-discrete model equations
The system of equations ( 1) is discretized using a nite volume method on a uniform staggered grid, sketched in Figure 2.This discretization naturally conserves mass for each uid separately, and momentum for both uids combined.The rst two components of (the phase masses) and the pressure are de ned at the centers of pressure volumes, which have a cell size of ∆ = ∕ .The last two components of (the phase momenta) are de ned at the centers of velocity volumes.
pressure cells velocity cells )∆ On this staggered grid we de ne a local discrete vector of unknowns as follows: The choice of using − 1∕2 in the de nition of instead of + 1∕2 is arbitrary.Note that 1, ( ) ≈ 1 ( , ) (and similar for the other entries); the notation is on purpose kept very close to the notation of the continuous model, but can be distinguished due to the extra index which the discrete variables carry.Another notable di erence is that the cell sizes are included in the discrete unknowns, so that they have units of mass and momentum.
The last equality in (34) describes the relations of the discrete conservative variables to the discrete primitive variables (cross-sections and velocities).Here we have introduced the following notation for interpolation operators [9]: The numerical scheme is implemented in terms of the conservative variables 1, through 4, −1∕2 , but the primitive variables can be extracted in post-processing according to the given relations.
The notation with as a discrete local vector of unknowns allows us to write the discrete scheme in vector form as Here, we have de ned −1∕2 as The numerical uxes and numerical pressure terms are left unde ned in this section, because we will de ne them based on the requirement of energy conservation, in section 5.
The pressure terms are non-conservative and are not written as the di erence between an in ow and an out ow of the nite volume cell.However, with the staggered grid employed here, one can see that 3, −1∕2 and 4, −1∕2 are directly and naturally connected to the pressure at the neighboring grid cells.Analogous to the incompressible (multi-dimensional) single-phase Navier-Stokes equations (for which staggered grids are known to lead to strong coupling), this pressure-velocity coupling is necessary to prevent checkerboard patterns, and would be much more di cult to achieve on a collocated grid.
Just as in the continuous case, these constraints are used to set up a Poisson equation for the pressure.The semi-discrete momentum equations are rst summed to obtain Expanding and substituting the de nition of −1∕2 yields After taking the di erence between this equation and the same equation for index + 1∕2, and applying (38), we obtain the discrete version of ( 11): System (36) is discretized in time using the fourth-order semi-explicit Runge-Kutta method described in [28].At each stage of the Runge-Kutta time step, a predictor-corrector algorithm is applied: the momentum equations are rst solved without including the pressure terms, the discrete Poisson equation is solved for the pressure using these intermediate momenta, and the momenta are updated in a projection step using the calculated pressure.This ensures that the volume and volumetric ow constraints are satis ed at all stages.We solve (41) iteratively, using a preconditioned conjugate gradient method.The time integration method is fourth-order accurate for all variables, and requires a restriction to the CFL-number based on the eigenvalues of the TFM.

Boundary conditions
In the case of periodic boundaries, the domain is divided into pressure volumes and = velocity volumes.There are no special boundary points: the scheme as laid out in subsection 4.1 applies everywhere, looping around the domain.
For closed boundaries, there are interior pressure points and = − 1 interior velocity points.The rst interior pressure node is located at = ∆ ∕2, the rst interior velocity node is located at = ∆ , and similarly for the last nodes at the end of the domain [28].For both the pressure and velocity grids, there are boundary points in addition to the interior points, one at each side of the domain.When calculating the discrete energy on the velocity grid (see subsection 5.1), it is important to include the half-volumes between the boundary points and the rst and last interior points.
At the boundary points, the mass uxes ( and ) are speci ed, and and follow via an analysis of the characteristics corresponding to the incoming and outgoing waves at the boundary.In the case of closed boundaries, as used in this work, the mass uxes are set to zero.Note that the characteristic analysis incorporates the volume constraint (37), and no boundary condition is needed for the pressure (the pressure at the boundaries has no in uence on the solution in the interior).For more details on the implementation of the boundary conditions we refer to [28].

Outline: conditions for discrete energy conservation
In the discrete case, just as in the continuous case, we want to satisfy a local and global energy equality.The use of a staggered grid instead of the commonly used collocated grid (e.g.[9,31]) makes it straightforward to obtain an energy-conserving discretization of the non-conservative pressure term, but introduces new challenges in terms of the de nition of the discrete local energy, which is not unique anymore.
We choose to de ne the local energy at the velocity grid points, i.e. we choose −1∕2 = ( −1 , ), and are aiming for a discrete version of (12): with ℎ = ℎ( −1 , ) and = ( ) as the numerical energy uxes.This choice means that the potential energy terms and 1 and 2 in the kinetic energy terms need to be interpolated, but 3 and 4 do not require interpolation.With this choice, we obtain energy-conserving expressions for 3, and 4, in a constructive manner (after choosing advantageous expressions for 1, −1∕2 and 2, −1∕2 ).It is also possible to de ne the energy at the pressure grid points, and obtain an energy-conserving discretization, but in that case it is necessary to substitute trial solutions for 3, and 4, , and interpolation of the pressure is required in the expression for (see the remark at the end of subsection 5.4).We would like to emphasize that (42) is not being solved as an additional equation; instead it will be shown to be a consequence of the discrete mass and momentum equations given in section 4, if the numerical uxes and are chosen appropriately.
If (42) holds, it can be summed over all nite volumes to yield where the last equality should hold in the case of periodic or closed boundaries.Here we have de ned the global discrete energy as the discrete counterpart of ( 14): Like in the continuous case, the art is to nd expressions for −1∕2 , ℎ and such that equation ( 42) is satis ed.In addition, the numerical ux −1∕2 needs to be constructed.We will outline the steps to obtain these quantities in a manner parallel to the continuous derivation in section 3.
First, we postulate an energy which will be based on the energy found for the continuous case.Note that the dependence could be expanded to additional grid points if required, but we will introduce an energy for which this is not necessary.Second, calculate the vectors of entropy variables, de ned as Here the rst index refers to the index of the energy, and the second index refers to the conservative variables to which derivatives are taken.
For the energy given by ( 45), the time derivative can be expressed as Here the brackets represent dot products over the vectors (at a certain grid point), just as in the continuous case.The right-hand side of equation ( 46) follows by substituting equation (36) for and − 1: where we have introduced the following notation for jump operators [9]: Comparing with (42) we see that the energy uxes ℎ and need to satisfy These conditions are analogous to ( 17) and ( 19) for the continuous case, with discrete jumps corresponding to derivatives with respect to .Together, conditions (49) and (50) guarantee that (47) can be written as (42), thus proving conservation of the discrete local energy (45).
The challenge is to nd the proper combination of discrete expressions for −1∕2 , ℎ , , and −1∕2 which are consistent approximations to their continuous counterparts in such a way that the local energy conservation equation is satis ed.This is a di cult problem, since we have multiple degrees of freedom ( −1∕2 , ℎ , , and −1∕2 ), and the solution might not be unique.To simplify the construction, we will use the concept of entropy potential introduced in subsection 3.3: after choosing a certain −1∕2 and −1∕2 , this yields straightforward conditions on the uxes −1∕2 to be energy-conserving.

Choice of discrete energy and energy uxes
In this section we propose an energy −1∕2 , and verify that this energy is conserved by the pressure terms of the discrete model (energy conservation for the ux terms is treated in subsection 5.3 and subsection 5.4).Recalling the continuous energy (21), we de ne a discrete energy −1∕2 = −1∕2 ( −1 , ): Other choices are possible because on a staggered grid interpolation is required, and the interpolation may be carried out in various di erent ways 3 .Our choice (51) is one of the most straightforward choices for the energy that is consistent with the continuous de nition, when the energy is de ned at the velocity grid points, and leads to an elegant form of the energy-conserving discretization (see also the remark at the end of subsection 5.4).
We use identities given in Appendix A to calculate the vectors.They are given by and and their sum is consistent with (22).The pressure terms in (47) need to satisfy condition (50), which can be rewritten to obtain the discrete version of (20): On a staggered grid, it is straightforward to satisfy this condition by choosing for since with this choice we have Consequently, condition (54) can be written with the volumetric ow constraint (38) as so that (54) (and (50)) is satis ed when is given by Note that our constraint-consistent time integration method enforces that the volumetric ow constraint is satis ed up to machine precision [28].

Reformulation in terms of the entropy potential and conditions on numerical uxes
The objective of nding energy-conserving numerical uxes is better served by reformulating condition (49) in terms of the entropy potential, because this results in an alternative, constructive, condition for nding energy-conserving uxes.The uxes are then based on the entropy potential −1∕2 instead of the energy ux ℎ .Similar to subsection 3.3, we rewrite the left-hand side of (49) as: which can be interpreted as a discrete version of the product rule ⟨ , ⟩ = ⟨ , ⟩ − ⟨ , ⟩.We have made use of the following de nitions: These de nitions are such that we only interpolate or take jumps between vectors with the same relative indices.Instead of directly choosing , it is more natural to use the last terms in (58) and de ne the jump in (similar to (28) since this leads to the following 'implied' de nition of : which is consistent with (27).The advantage of (61) over ( 49) is that we have a condition on the ux itself, rather than on the jump in the ux.Once −1∕2 and −1∕2 have been chosen and −1∕2 has been derived, ℎ can be determined from ℎ = ⟨ , −1∕2 , −1∕2 ⟩ + ⟨ , +1∕2 , +1∕2 ⟩ − .
We propose now the following discrete entropy potential for the equations: This is a straightforward discretization of (29).Given the expressions for ((52) and (53)), condition (61) can now be evaluated to yield the numerical uxes .In order to be able to derive from the (scalar) condition (61) multiple equations for the individual numerical uxes, we split 3, −1 and 4, −1 into an advective component (denoted by subscript ) and a level gradient (or gravity) component (denoted by subscript ): As a consequence, condition (61) can be split into the following four separate conditions by collecting terms featuring and those not featuring , and by using the functional dependencies assumed for the uxes: 3, , = 0, (65a) , .
These conditions have been obtained analogously to their continuous equivalents (31).In the continuous case, the uxes were known and these conditions were satis ed by construction.In the discrete case, these conditions will be used in the next section to determine the numerical uxes.

Derivation of energy-conserving numerical uxes for the TFM
System (65) is a system of four equations for six unknowns.To nd a solution we assume, based on the continuous expression, that This choice is motivated by the fact that it requires no interpolation, and moreover is such that the discrete Poisson equation (41) follows naturally from the discrete volumetric ow constraint (38) (which is used in our time integration method [28]).Substituting 1, −1∕2 in (65a) and 2, −1∕2 in (65b) yields directly 3, , = To get the gravity component of 3, and 4, , substitution of 2, −1∕2 in (65d) leads to .
After signi cant rewriting, this yields the following expression for the gravity component of 4, : and a similar expression holds for 3, , .
The rst term on the left-hand side is easily recognized as the discrete counterpart of − .In order for the discrete expression to be practical and match the continuous expression, the second term must vanish, and we require the following conditions to be satis ed: In the continuous case a continuous version of these conditions, given by ( 26), is also required, and these can be shown to be satis ed exactly via manipulation of the continuous derivatives.The same manipulation is not possible with discrete jumps, so that in the discrete case these conditions are not satis ed in general, and the second term in (68) does not generally vanish.This means that we cannot obtain a practical energyconserving discretization for arbitrary geometries (at least not with the conventional staggered-grid nite volume method that we have employed).Even though conditions (69) are not generally exactly satis ed in the discrete case, we can show that they are approximately satis ed for arbitrary duct geometries, and that they are exactly satis ed for speci c geometries such as a channel.This can be shown by evaluating both sides of (69) using Taylor series.We expand , −1 and , −1 into Taylor series around = , , and expand , −1 around = .These Taylor series are combined to obtain expressions for , −1∕2 , , −1∕2 , and , −1∕2 .With these expressions the left-hand side of (69) evaluates to where (.) indicates (.) evaluated at , .The right-hand side of (69) evaluates to At this point we apply relation (A.6) from Appendix A to the discrete quantities used here: and from this we can derive d 2 , and Substitution of these relations in (70), and comparison of the result to (71) yields This derivation can be carried out with similar results for the upper uid.These relations show that for arbitrary duct geometries, the geometric conditions (69) are satis ed only approximately in the discrete case.This stands in contrast to the continuous case, where the equivalent geometric conditions are satis ed exactly (for arbitrary geometries).
Fortunately, for a 2D channel geometry d ∕d = 1 and d 2 ∕d 2 = 0, and all higher order derivatives are zero, so in this case (69) is exactly satis ed.This means that the 2D channel geometry is an important special case for which we obtain the following numerical uxes: These uxes are energy-conserving for other geometries with d 2 ∕d 2 = 0, but not for geometries with curved sides, such as the pipe geometry.
The nal collection of energy-conserving numerical uxes is given by ( 66), (67), and (73).Of these, (66) and ( 73) are locally exact, and (67) involves second order accurate central interpolation.Together they form the numerical ux vector Here the ux is rendered in terms of primitive variables only for ease of interpretation; the implementation of the numerical ux (and of the discrete energy) in the numerical code is completely in terms of the conservative variables.
Remark 1.The di culty to satisfy condition (69) for arbitrary cross-sectional geometries is not dependent on the choice of −1∕2 , nor is it due to the interpolation of the potential energy to the velocity grid points (as needed on a staggered grid).This is shown in Appendix B by applying a global energy analysis.
Remark 2. The proposed discrete energy ( 51) is a consistent approximation to (21) which is conserved by the numerical uxes given by ( 74).However, it is not unique.For example, an alternative de nition is In this formulation the energy is de ned on the pressure grid, and the energy conservation conditions and local energy conservation equation can be adapted to accommodate for this.With a similar change in the entropy potential, it is again possible to derive a set of energy-conserving numerical uxes, which turn out to be the same as those given by (74).As the issue of the geometric relations also persists with this choice, there seems no clear advantage over our proposed formulation.

Numerical experiments
We perform numerical experiments for a 2D channel geometry, with the goal of verifying conservation of the discrete global energy, as discussed in subsection 5.1: The model for which we perform the experiments will not include source terms such as wall friction and interface friction, or di usion, since these would lead to dissipation of energy in the continuous analysis.The test cases are chosen such that no discontinuities appear, for which the continuous analysis is invalid, since this would also necessitate dissipation of energy.Furthermore, the numerical experiments performed in this section will all be in the 'well-posed regime' of the TFM, meaning that the initial conditions are chosen such that the eigenvalues of the model are real, and remain so.We use the discretization as outlined in section 4, with the numerical uxes given by (74).The vector of the pressure term is given by (55).We noted earlier that the scheme is spatially exactly energy-conserving, but not temporally.However, we can still obtain energy conservation by taking the time step su ciently small.The di erence between the initial energy 0 ℎ and the nal energy ℎ after time steps should then be in the order of the machine precision, and we shall term this di erence the 'energy error'.

Gaussian perturbation in a periodic domain
We consider a test case with periodic boundaries, so that e ectively we do not need to take the boundaries into account.We introduce a perturbation in the hold-up = ∕ of the form , with ∆ = 0.2 and = ∕10, and the length of the domain.This produces a Gaussian perturbation centered at the middle of the domain.The initial velocities are left at zero, which ensures exact initial satisfaction of the volumetric ow constraint (8) (in fact, = 0).We use parameters similar to those used in the Thorpe experiment [32], as described by [21].They are given by Table 1.The choice for a large upper uid density is deliberate: it ensures that all terms in the expression for , (51), are signi cant.We employ = = 40 nite volumes with ∆ = ∕ and let the simulations run until = 30 s, with ∆ = 0.001 s.The perturbation splits symmetrically into a left-traveling and a right-traveling wave, which travel through the periodic boundaries, to eventually come together in the middle and reform the initial perturbation approximately.We show the evolution of the hold-up and velocity in Figure 3, roughly up to the point that the waves meet at the boundaries of the domain.In this test case we have a signi cant exchange between kinetic and potential energy, which can be seen in Figure 4 (left panel).The total energy is conserved up to machine precision, as can be seen in the right panel of the gure.The mass of each phase and total momentum are also conserved, and the volume constraint and volumetric ow constraint are satis ed, up to machine precision (see also [28]).As time progresses, nonlinear e ects start to play a role, leading to more irregular behaviour of the potential and kinetic energy as a function of time.The sum of the two stays exactly constant, con rming our theoretical derivations, and showing that our newly proposed numerical uxes for the TFM lead indeed to an energy-conserving discretization method.
We give further evidence that the energy is conserved exactly by the spatial discretization, and limited only by a temporal error, by plotting the convergence of the energy error with re nement of the time step.Figure 5 shows a fourth order convergence rate with ∆ , in agreement with the fourth order accuracy of the Runge-Kutta time integration method.The convergence continues up to machine precision, which is reached around ∆ = 0.001 s, as was used for the results in Figure 4, con rming that the spatial discretization conserves energy up to machine precision.

Sloshing in a closed tank
We now consider a test case with closed (solid-wall) boundaries, for which energy conservation is expected to hold because the uxes ℎ and involve multiplication with 3 and 4 , which are zero at the boundaries.The test case features a closed rectangular tank in which the two uids are brought out of equilibrium, so that sloshing occurs.The parameters are identical to those of the previous test case, see Table 1, except that the initial condition for the hold-up perturbation is di erent.It is given by with ∆ = 0.2.This yields a straight slanted interface, with = 0.3 at the left boundary and = 0.7 at the right boundary: see Figure 6.This is not a typical sloshing case, since the TFM was designed to model long-wavelength phenomena, and indeed we have taken ≫ .Therefore we are not able to explicitly capture typical sloshing phenomena such as wave breaking.However, the e ect of such small-scale phenomena on the averaged ow may be included in the model via closure terms [13].With accurate closure terms, the TFM can closely match DNS results, as shown in [5].This is not included here, as this would lead to dissipation of energy and not allow us to show the energy-conserving properties of our proposed numerical discretization.
Like in the rst test case, initially the total energy of the system consists of only potential energy.Under the presence of gravity (via the level gradient terms) the interface starts to atten, which is achieved via a right-running and a left-running wave, that emanate from the left and right boundary, respectively.Around = 7 s the interface is almost completely at, and all potential energy has been converted into kinetic energy, and the interface starts to slant ('slosh') again in the opposite direction.Figure 6 shows this behavior up to approximately the point that the lower uid reaches its maximum height at the left boundary.Note that the evolution of the hold-up fraction is not exactly symmetric, amongst others because the wave speed in the 'deep' part is di erent from the wave speed in the 'shallow' part.Also in this test case, the mass of each phase is conserved up to machine precision, but there is a (physical) in ow of momentum at the boundaries, due to the level gradient terms.Figure 7 shows the exchange of potential and kinetic energy as a function of time.Similar to the previous test case, exact energy conservation is achieved with our proposed spatial discretization, if the time step is ne enough (here ∆ = 0.005 s, and = 40).If the time step is not ne enough, a (small) energy error is made, which converges with fourth order upon time step re nement, as is shown in Figure 8.The ability to conserve energy in this closed system is an important step in order to obtain delity in the simulation results.Non-energy-conserving schemes, e.g.schemes that dissipate energy, would introduce arti cial (numerical) damping of the sloshing movement and incur a loss in the liquid height reached at the boundaries.In a way, the sloshing movement can be compared to a moving pendulum [25], for which it is well-known that conservation of the total energy (the Hamiltonian) is an important property that should be mimicked upon discretization in order to achieve realistic long-time behavior.

Traveling wave
Finally, we perform a test case with a traveling wave in a periodic domain.The ow is uni-directional and strati ed, with a velocity and density di erence between the two uids.We consider a steady base state, upon which a small periodic perturbation is introduced, of which we study the evolution in time.This case is similar to test cases examining the Kelvin-Helmholtz instability, such as in [19,28].However, here the perturbation will be stable since the ow is inviscid and in the (linearly) well-posed regime.
Most of the parameters are again identical to those given by Table 1, but the initial conditions for the hold-up and the uid velocities are di erent.We set ,0 = 0.4 and ,0 = 1.For we take ,0 = 1.187 (which is the value that would result in a steady ow with wall and interface friction 4 ).
In order to construct an initial perturbation that results in a traveling wave, we conduct a linear stability analysis of the TFM [19].The analysis is conducted in terms of its primitive variables in the form As exact solutions we obtain waves of the form with ∆ the amplitude of the perturbation in each variable.The relative amplitudes in ∆ are such that ∆ is an eigenvector corresponding to one of two dispersion relations ( ).
The initial perturbation is de ned as (77), with = 0. We take a wavenumber of = 2 ∕ m −1 and calculate the corresponding angular frequencies, of which one is selected.The chosen mode is = 3.982 s −1 .
Setting ∆ = 1 ⋅ 10 −2 , the amplitudes of the other variables are calculated so that ∆ is an eigenvector corresponding to this mode: (∆ ) = 1.00 ⋅ 10 −2 3.99 ⋅ 10 −3 4.51 ⋅ 10 −3 −2.30 .This ensures that the other mode is not present in the initial perturbation, so that we can study the isolated behavior of one mode.A projection step is then performed in order to make the initial condition satisfy the constraints (see section 4).
The initial condition is shown in Figure 9, along with its evolution in time, which is computed up to = 30 s. Setting the initial condition this way yields a wave traveling to the right at velocity ∕ = 1.16 m s −1 , which remains of approximately constant amplitude since the ow is inviscid and in the well-posed regime, so that has no imaginary component.The traveling wave can deform due to the nonlinear character of the governing equations, which is neglected in the linear stability analysis.This is made apparent by the snapshots of the solution shown in Figure 9, which are separated by an integer number of wave periods: at the time of the last snapshot the wave has traveled through the domain 18 times.The solutions do not completely overlap and we see wave steepening taking place.
Figure 10 shows the evolution of the energy.In this case, the exchange between kinetic and potential energy is small relative to the total energy of the base state.This is due to the fact that the wave is roughly constant in time, up to a displacement which does not change the energy.
The total energy can again be seen to remain constant up to a high precision.Like before, this is achieved by using a small time step (∆ = 0.005 s), with a modest spatial resolution ( = = 40).Figure 11 shows how the energy converges with time step re nement.The convergence rate is fourth order over a wide range of time steps (matching the order of the time integration method), demonstrating that also for this test case, the spatial discretization conserves energy.While the solution moves away from the stable traveling wave predicted by linear analysis, its energy remains constant with time.

Conclusions
In this article, we have derived the result that the total mechanical energy (sum of kinetic and potential energy) is a secondary conserved quantity of the incompressible and isothermal TFM.This result is in line with the well-known fact that multi-dimensional incompressible frictionless ow equations conserve mechanical energy.Our novel insight is that this conservation statement still holds after averaging: the averaging procedure used to obtain the 1D TFM does not interfere with the energy conservation property.The approach was based on the formulation of entropy variables and an entropy potential, similar to what is commonly done for the SWE, but with two main di erences: (i) we have included a non-conservative pressure term in our analysis, which is shown to be energy-conserving, and (ii) we have obtained our results independent of the duct geometry, which may be a 2D channel or a circular pipe, or any other closed cross-sectional duct shape.
The second novel result of this paper is a set of numerical uxes that conserve a discrete form of the mechanical energy.A discretization on a staggered grid was proposed in order to keep the energy conservation property of the non-conservative pressure terms in a discrete sense.Although the use of a staggered grid implies that the choice of a discrete energy and entropy potential is not unique, we were able to propose a combination which is such that the discrete analysis is consistent with and analogous to the continuous analysis.However, one important di erence between the continuous and discrete cases remains, namely in the analysis of the level gradient terms (for arbitrary geometries).A geometric relation between the potential energy and the interface height is satis ed exactly in the continuous case, but only approximately in the discrete case.Fortunately, for the speci c case of the 2D channel geometry, the condition is satis ed exactly, and the discrete level gradient reduces to a form which parallels the continuous form perfectly.For other geometries, such as the pipe, a small numerical energy error persists in the discrete analysis.
Our theoretical derivations are supported by numerical experiments, which show that the proposed energy is indeed exactly conserved by our new spatial discretization in both periodic and closed domains.Building on previous work [28], the discretization also conserves mass and momentum, has strong coupling between momentum and pressure, and is constraint-consistent.In these experiments the temporal error was negligible (due to a combination of high-order time integration and small time steps), but for future work it is suggested to also make the time integration method energy-conserving [26].Furthermore, the e ects of wall friction and pipe inclination need to be added into our formulation.
Our energy-conserving formulation of the TFM provides a foundation for investigating the nonlinear stability of the model.For related models, the energy acts as a norm or a convex entropy function of the solution, providing stability bounds, and it should be investigated if the TFM energy has similar implications.While we have only considered smooth solutions, in the theory of entropy stability [30], energy is dissipated at discontinuities.Therefore, in order to deal with shocks it seems necessary to add suitable di usion to our formulation, so that the energy becomes strictly decreasing [6].
The integrals (5) which appear in the governing equations of the two-uid model are geometry-dependent: with (ℎ) the local width.Note that ( ) = int .For a 2D channel geometry, with = and = , the width is given by (ℎ) = 1 and the integrals evaluate to where we have substituted = − .For the pipe geometry, we make the transformation ℎ = (1 − cos ( * )), with * the integration variable and the wetted angle, to get [27] The following derivatives of and are needed in order to calculate ∕ : We use Leibniz' rule to calculate which can be evaluated by substituting the expressions for and .In order to calculate as given by ( 22), we need the derivatives d ∕ and d ∕ .They are found by di erentiating (A.8) and (A.9), yielding

Appendix B Global energy analysis
The main text has described a way to derive the local semi-discrete energy conservation equation given by (42).In the case of periodic or closed boundaries, this can be integrated in space to yield global energy conservation.In this section, we directly derive the global energy conservation equation without the intermediate step of the local energy.This allows us to skip the step of choosing an entropy potential, which means that the derivation will contain less assumptions.On the other hand, the obtained conditions on the numerical uxes are not constructive, because they are conditions for the 'jumps' of the numerical uxes, rather than for a single numerical ux at one discrete point.
The scheme (36) described in section 4 for a certain pressure volume and velocity volume − 1∕2 can be extended to describe the evolution of the entire state vector ℎ : We split this condition into two conditions: one proportional to and one not proportional to : ⟨ ℎ , ℎ ⟩ = ⟨ ℎ , ℎ ⟩ + ⟨ ℎ , ℎ ⟩ .
The advective condition is given by Here, the sum over the entries on the rst lines evaluates to zero, since each term has a matching term of opposite sign and index shifted by 1 (even the boundary terms, in case of periodic boundaries).In order for this to also hold for the terms in the second and third lines, we need to satisfy the condition Now, in order for this to be conservative, we need the rst term in each line to be equal but opposite in sign to the second term in each line (shifted in index by 1).This yields the following conditions: which upon substitution of (73) reduce to the geometric conditions (69).
In conclusion, the results of the global discrete analysis are consistent with our local discrete analysis.The additional insight from the global analysis is that the geometric conditions (B.6) or (69) are independent of the choice of the entropy potential.This con rms that the choice of entropy potential does not limit the results.

Remark 3.
The global energy analysis can also be performed without requiring interpolation of the potential energy to the velocity grid points, as needed in the de nition of −1∕2 given by (51).Instead, one can directly de ne It can be veri ed that this leads to the same ℎ as given by (B.5), and consequently the geometric condition (B.6) remains present.

Figure 1 :
Figure 1: A schematic of strati ed two-uid ow in ducts (a circular pipe segment is shown as an example) described by the 1D TFM.

Figure 3 :
Figure 3: The initial evolution of the Gaussian perturbation, up to the point that the boundaries are met.Left: lower uid hold-up.Right: lower uid velocity.

Figure 4 :
Figure 4: Conserved quantities for the Gaussian perturbation test case.Left: potential, kinetic and total energy relative to their initial values.Right: ℎ − 0 ℎ ∕ 0 ℎ .

Figure 5 :
Figure 5: Convergence of the energy error ℎ − 0 ℎ ∕ 0 ℎ with time step, for the Gaussian perturbation test case.

Figure 6 :
Figure 6: The initial evolution of the sloshing simulation, approximately up to the point that the lower uid reaches its maximum height at the left boundary.Left: lower uid hold-up.Right: lower uid velocity.

Figure 9 : 13 Figure 10 :Figure 11 :
Figure 9: The initial perturbation travels to the right with time.Consecutive snapshots are separated by a time interval of 6 wave periods.Left: lower uid hold-up.Right: lower uid velocity.
of the derivatives appearing on the right-hand sides can also be evaluated using Leibniz' rule: geometric quantities are used in(21) and de ned as:

Table 1 :
Parameters for the Gaussian perturbation test case.