On the discrete equation model for compressible multiphase fluid flows

The modeling of multi-phase flow is very challenging, given the range of scales as well as the diversity of flow regimes that one encounters in this context. We revisit the discrete equation method (DEM) for two-phase flow in the absence of heat conduction and mass transfer. We analyze the resulting probability coefficients and prove their local convexity, rigorously establishing that our version of DEM can model different flow regimes ranging from the disperse to stratified (or separated) flow. Moreover, we reformulate the underlying mesoscopic model in terms of an one-parameter family of PDEs that interpolates between different flow regimes. We also propose two sets of procedures to enforce relaxation to equilibrium. We perform several numerical tests to show the flexibility of the proposed formulation, as well as to interpret different model components. The one-parameter family of PDEs provides an unified framework for modeling mean quantities for a multiphase flow, while at the same time identifying two key parameters that model the inherent uncertainty in terms of the underlying microstructure.


Introduction
The dynamical evolution of two (or more) distinct phases (of matter) is often referred to as multiphase flow and it is a very important topic of study in a broad variety of engineering systems, even though it is by no means limited to modern industrial design and can be observed in many natural/biological phenomena. A very limited list of references for multiphase flow include [17,21,15,44,6,7,36,5] and references therein.
The simplest, yet very representative, form of multiphase flow is two-phase flow. The mathematical modeling of two-phase flow arguably originated in the so-called multi-fluid models. Herein, one assumes that the dynamics of compressible inviscid fluid mixtures is modelled by the Euler equations [36], where the characteristic middle field (contact discontinuity) consists of a material interface, if the adjacent data belong to different phases. Different parameters in the equations of state (EOS) are introduced in these models to represent the inherent heterogeneities in terms of the discontinuous variation of the pressure-density relations. Finally, additional conservation laws are inlucded to model species advection [37,29,1,2,3,45,25,11]. Despite the inherent simplicity and flexibility of this approach, such models are often marred by spurious velocity and pressure oscillations near material interfaces [1,25,2,3], excessive numerical diffusion [43], when approximated via classical schemes and negative mass fractions [29].
An alternative and more popular approach, based on the theory of multiphase flows [21,15], assumes each phase to be distinct and described by its own set of equations, typically the Euler equations. Pioneering works in this direction include those of Stewart and Wendroff [44] and Bear and Nunziato [6], see also [36]. This approach has now been extended into a wide variety of possible models. Following the observation that different phases interacts through the interface up to reaching uniform conditions [7] (i.e. they move with approximately the same pressure and velocity), the resulting set of equations is classified according to the set of independent variable they consider [49]. Restricting the discussion here to one space dimension, we start with the so called four equation models [24,48], which essentially resemble the reactive Euler equations and lead to similar difficulties as those experienced with the multifluid approach described above.
Next, one considers the so-called five equations models [7,24,33,26,42], where one assumes a fully mechanical equilibrium between the phases,implying that the mixture is macroscopically moving with one-pressure and one-velocity. In [33], it is shown how to derive the five equation models from the Baer and Nunziato one by a formal asymptotic expansion assuming that the relaxation parameters tends together toward infinity, while their ratio stay bounded. In the case of non smooth solutions, a set of jump relations for the five equation model was provided in [42].
One can follow [36] and relax the assumption of mechanical equilibrium across phases. The resulting seven equation model requires the introduction of stiff source terms to model the underlying thermodynamics and leads to the removal of spurious oscillations around material discontinuities. Moreover, the source terms force a relaxation to a single pressure and velocity recovering an experimentally observed fact in two-phase flows. Moreover, the zero relaxation limit of these models results in the five-equation model of Kapila et al. [24].
Inspite of the tremendous progress made with regards to the modeling of two-phase flows as described above, several pressing issues remain. To start with, these mathematical models involve nonconservative products which make conservation of energy potentially difficult. Moreover, a mathematically sound solution concept, together with rigorous proofs of well-posedness, even in one space dimension, is extremely challenging. Notable exceptions are presented in [23,27,34] where the authors provide a rigorous mathematical treatment of a simplified version of the Baer-Nunziato equations.
Furthermore, from a modelling perspective, a stark shortcoming of the many of the afore-mentioned models lies in the fact that the interfacial velocity and pressure are difficult to determine, see [20,32,8,36,30,13,35,6,12,10,14,40,19,41] and references therein for a discussion of this issue as well as possible remedies.
Given these shortcomings of the afore-mentioned models, one can see that there is no consensus on what constitutes a suitable modelling framework for two-phase flows. In particular, an uniform description of the vast range of flow regimes, ranging from isolated interfaces to fogs and microbubbles, within the purview of a single predictive model is extremely challenging. The search for such a framework brings us to the so-called Discrete Equation Method (DEM) of [5], see also [2]. Inspired by the Godunov method and well-established theories of ensemble averaging [15], DEM entails the statistical description of each phase in terms of its own equation of state and allows for, in principle, all possible flow regimes. A multiscale formulation allows one to incorporate information from finer scales. One can think of DEM as a mesoscopic model as its does not require an explicit description of the underlying microstructure.
Despite its promise as a suitable modeling framework for multiphase flows, DEM still requires userdefined ansatz (closure relations) on the probability coefficients that arise in course of the ensemble averaging procedure. Although many papers such as [38] suggest modifications for overcome this issue, for instance in the case of simulating dense-to-dilute transitions by coupling the underlying Euler equations with an evolution equation on the number of dispersed particles, it is fair to say the design of a flexible general purpose DEM type model, which can describe various flow regimes is still outstanding.
These limitations of the DEM approach constitute the starting point of the current paper. Herein, we will carefully develop and analyze the DEM approach for describing two-phase flows in one space dimension, while neglecting heat and mass transfer. Our main aim would be to characterize the probability coefficients that arise in the DEM framework of [5] such that all possible flow regimes can be described by DEM. This will allow us to encapsulate all phase interactions in terms of a single parameter that interpolates between disperse and stratified flows. Moreover, simple relaxation procedures will also be investigated. This will allow us to study numerically, how different choices of parameters leads to the recovery of different flow regimes, enabling a thorough analysis of the expressivity as well as limitations of DEM for different regimes of multiphase flow.
The rest of the paper is organized as follows: In Section 2 we summarize the DEM procedure, highlighting the modelling assumptions related to such procedure. Section 3 is dedicated to the analysis of the probability coefficients resulting from the previous section, and Section 4 derives the corresponding one-parameter limit along with the numerical strategy to solve it. Finally, Section 5 include the numerical experiments we have performed on such models, and discussion of the outcomes is carried out in Section 6.

The Discrete Equation Method
In this section, we will present the discrete equation method for modeling two-phase flows in one space dimension. We start with a succinct presentation of the ensemble averaging theory on which DEM is based.

The ensemble averaging theory
In the following we recall the procedure of [5] for a biphasic Eulerian flow without mixing. Phase transition is excluded from the present study and we suppose that heat transfer is too slow compared to mechanical relaxation [24]. We consider two phases Σ 1 and Σ 2 , each governed by the Euler equations The notation is classical: ρ (k) , u (k) , p (k) denote the density, velocity and pressure of the phase k ∈ {1, 2}. The total energy E (k) = 1 2 u (k) 2 + e (k) , where e (k) denotes the internal energy. Different choices of equation of state (EOS) have severe implications on the flow regime and a typical issues in multiphase flow is the determination of a methodology that handles different EOS.
As it is well-known [15], a prime characteristic of multiphase mixtures is that there is uncertainty in the exact location of the particular constituents at any particular time. In turn, from the practical point of view, for a given set of initial and boundary conditions, a single measurement of such experiment carries limited information about the mean and distribution of dispersed particles that generated such results. For this reason, modern multiphase flow theory is described in averaged sense. In our case, we aim at considering both the spatial rearrangement of disperse particles and the statistical description of repeated sampling for a fixed set of initial and boundary condition.

Notation
We hereby introduce some notations. Let (Ω, F, P) be a probability space on R. We denote the physical space of interest (i.e. domain) by an open set D ⊆ R d , where d ∈ N is the spatial dimension. The time horizon is denoted by T > 0, and the any time considered for our simulations is denoted by t ∈ [0, T ]. We aim at including the randomized dependency of quantities of interest by taking random fields between the spaces (Ω, F, P) and the space of p-integrable functions L p (D × R + ; U ), with U ⊂ R N . Here N ∈ N is the number of quantities of interest of the system under consideration.
Existence and uniqueness (well-posedness) of solutions for systems of hyperbolic conservation laws is restricted to one-dimensional (d = 1) and for sufficiently small initial data [9]. More sophisticated solution paradigma [18] are also available but they are out of the scope of this work. We therefore restrict our description to the cases d = 1.
In such a case, weak -solution are typically found in the subspace BV (D × R + ; U ). We will consider random variables between the spaces (Ω, F) and (X , B(X )), where the topological space X = L 1 is endowed with the Borel-sigma algebra B(X ), as to make each continuous function measurable. Let ω ∈ Ω be a fixed realization. At each time level t ∈ [0, T ] we will assume that there exist a pair of open sets D 1 (t; ω), D 2 (t; ω) affected by only one phase, namely D k (t; ω) := {x ∈ D | phase k is present at (x, t)} such that 1. (Non-mixing condition) Only one phase is present at each space-time location: 2. (Saturation condition) No vacuum is generated at any space-time location: where ∂D denotes the frontier of D.
The interface between the two-phases is then defined according to the following relation: We introduce the characteristic function X (k) : Ω → X associated to phase k as the indicator function over the points of the domain D affected by phase k, namely Using standard theory of distribution, the characteristic function can be shown to satisfy the following topological equation (suppressing ω-dependence for notational convenience) [15] ∂ t X (k) + σ∂ x X (k) = 0 where σ is the interface velocity of the realization highlighted by X (k) (·; ω). Hence, one can also show that upon multiplication of (1) by the characteristic function it holds where the Lagrangian flux F (k) lag := F (k) and the subindex I denotes the interfacial value from the k-th side. We introduce the ensemble average operator E [15] that is assumed to commute with time and space derivative operators (these are commonly referred as Gauss and Leibniz Rules, which hold for well-behaved input functions). Taking ensemble average on (4) and (3), one obtains the following equation We thus introduce the notation that will be used throughout this paper: let where the ensemble average quantities are defined via so that E k := 1 2 u 2 k + e k . Using this notation the ensemble-average flux can be written as where F 0 k denotes the kinetic fluctuation of momentum and energy, that will be neglected in the following.

The DEM for Eulerian biphasic flow
Using the notation introduced in the previous section, the DEM method applies to the discrete setting: we consider a computational mesh (x i ) i=1,...,M ⊂ R and the associated control volume C i = According to the definition of Lagrangian Fluxes, one needs to identify/be able to compute the speed of the interface separating different components. This translates at the numerical level to the necessity of considering Riemann Solvers able to compute a contact-discontinuity σ. Given two initial states U L , U R we assume the solution of a Riemann Problem with possibly different phases at each side of the discontinuity to generate three waves (shocks or rarefactions separated by a contact discontinuity), in complete analogy to the single-phase theory. Given U L , U R ∈ R m , the speed of the contact-discontinuity/material interface is denoted by σ LR := σ(U L , U R ), while F (U L , U R ) and U (U L , U R ) denote the numerical flux and the numerical solution generated by solving the Riemann Problem with initial states U L , U R . The concrete forms of the numerical operators F, U depend on the Riemann Solver under consideration, for which popular choices are the HLLC or the Roe Riemann Solvers [46,. At each time step the preliminary stages of the method proceed as follows: at each time level t = t n , we have 1. Subdivide randomly the computational cell where ω aims at indexing the specific realization of X (k) .
2. Assign randomly in each subcell [ξ j , ξ j+1 ] the phases Σ 1 or Σ 2 with the state U (1) or U (2) . Up to merging adjacent subcells affected by the same phase, we have that within a volume two adjacent subcells contain different phases. We denote the interface velocity originating at the subnode ξ j as σ j , see Fig. 1.
Notice that the evolution of phase k ∈ {1, 2} Figure 1: Schematic representation of the prototypical generation of interfaces in the control volume can be written as where the Lagrangian fluxes F lag j := F (k) I j are evaluated at the only side affected by phase k of the interface moving with velocity σ j .
3. Obtain a semi-discrete approximation of the realization according to a Godunov type scheme: we approximate the flux integrals and the Lagrangian flux integrals by means of a Godunov type scheme for any r ∈ [t n , t n + s]. The notation X (k) j stands for the jump across the j-th interface moving with velocity σ j . Under the above assumptions and upon division by s in (10), the scheme reads Due to the alternate character of the distribution of data in the interior of the volume C i , one obtains the following relations: let us define the number of interior interfaces N int = N − 1 ≥ 0, then (a) N int is even : one can rearrange the summation as to arrive to (see Table 1) where the characteristic function χ over the even {Y = 1} is defined as Hence, by putting together the two instances that may occur, one ends up with where the perturbation variable θ (k) is defined as We thus assume that E θ (k) = 0 for each k, thus implying that the perturbation with respect to the first term N int i ) generated by an odd number of internal contributions is negligible in mean. Such an assumption is clearly not verified for a low number So the semi discrete scheme reads 4. Ensemble average of all realizations: taking ensemble average in (14) and with reference to the notation (6), we obtain, 3 The one-parameter mesoscopic scheme In order to be of practical use, the scheme (15) requires the specification of four different terms: : the average number of internal components of the dispersed phase in cell C i ; : the conservative numerical flux; 0 the left non-conservative term; : the right non-conservative term.
Building upon the work of Abgrall and Saurel [5], the aforementioned ensemble averages can be simplified by noticing that the random variable X (k) is in fact discrete, and its average can be written as the sum of all the instances multiplied by their probability of occurrence. In the following we will make use of the following notation for each phase index p = q ∈ {1, 2}, with the notation a prescribed time level t = t n . Notice that, these probabilities are defined in terms of different characteristic functions X (p) . Nevertheless, fixing the phase k = l ∈ {1, 2}, one can equivalently rewrite these latter probabilities in terms of one characteristic function Moreover, we define the flux indicator function and the notation a + := max(a, 0), a − := min(a, 0). Estimation of three of the above quantities is accomplished as follows: • Conservative Terms: We require that the Godunov state U * i+ 1 2 (0) [46] (i.e. the solution of the Riemann Problem at the right cell interface at time t = 0) belongs to the phase k or not -see Fig. 2. Hence, • Right Non-Conservative Terms: We need to make sure that a Lagrangian flux exists at the right interface, i.e. inflow is occurring. Figure 2: Schematic representation of the Godunov state U i+ 1 2 (t)| t=0 at the right cell interface.
• Left Non-Conservative Terms: We need to make sure that a Lagrangian flux exists at the left interface, i.e. inflow is occurring.
Using formulas (18)- (20), the final scheme reads Notice that the topological equation for the volume fraction is then recovered from (21) by formally choosing F ≡ 0 and U k ≡ 1, so that the Lagrangian flux reduces to F lag = 0 − σ · 1 = −σ. Hence, the numerical scheme is then of practical use, once the probability coefficients are defined.
Originally such probability coefficients were given by means of an ansatz, leading to a limited model, even if thermodynamically consistent [38,43]. We are going to close the model by proving convexity of such probability coefficients. The following proposition summarizes the properties each probability coefficient has to verify [5].
be a pair of probabilities coefficients defined in (16). Assume that Then the following consistency conditions must hold: for each p = q ∈ {1, 2} where the two neighbouring volume fractions verify the saturation condition A probability pair P will be termed a consistent probability pair if it verify (22)-(23), under the assumption that (24) holds.
Remark 1. Notice that the pair consisting of the two bounds in (23) is itself a consistent probability pair.
Due to this remark, Abgrall and Saurel proposed the following approximation for the probability coefficients An interesting features of this choice is that it has both mathematical and physical implications. First, from the mathematical point of view, it can be shown that fixing the probability coefficients then, as to form a consistent probability pair, we have no other choice The viceversa also holds. Furthermore, the pair P 0 constitutes an upper-lower bound for any pair of probability coefficients P i+ 1 2 , respectively, according to (22a) and (22b). Thus, such probability pair is an extreme point in the space of consistent probability pairs.
On the other hand, there is an interesting example that helps understanding the physical implication of choosing P i+ 1 2 = P 0 is not zero, that is, it is not zero the coefficient of each flux of the type F (Σ m , Σ m ) and F (Σ m , Σ n ) with m = n ∈ {1, 2} in (15). Figure 3: Schematic representation of a stratified flow at the numerical level.
By computing the corresponding probability pairs, one can convince oneself that these acts as fluxweights in (15), corresponding to the area of the surface through which a specific flux is applied. Unfortunately, the same computations would also be carried out in the case of disconnected phases at the interface, see Fig. 4. This second case will be called dispersed flow, where the phase that does not share a segment of the cell interface is called the dispersed phase. In such a case, each probability coefficient would still be non vanishing, thus introducing in the computation a non-zero numerical flux for the disperse phase, even if the regime is discontinuous. At the physical level, for the dispersed phase, this is equivalent to saying that a sound wave propagated in the dispersed phase in cell i gets propagated into the corresponding phase of cell i + 1 even if no material connection between the phases exists. Figure 4: Schematic representation of a numerical dispersed flow of phase k into phase l: the DEM method These considerations, motivated us to investigate the structure of such probability coefficients: the following proposition identify the complementary lower-upper bounds.
be a pair of probability coefficients defined in (16). Then for every p = q ∈ {1, 2} it holds under the assumption that (22) holds. Moreover, the probably pair is a consistent probability pair, provided that the saturation condition (24) holds.
Proof. See Appendix A.
Remark 2. Following the considerations made before Prop. 2, one can see that the probability pair P 1 is associated with a dispersed flow regime, see Fig. 4.
It light of this final remark, it becomes not surprising the following theorem.
Proof. See Appendix A.
Previous result leads us to substitute the probabilities appearing in (21) with the right hand side of (27), leading to a one-parameter semi-discrete scheme, modeling the mesoscopic description of the underlying two-phase flow.

Continuous limit and solution strategy
By suppressing the dependency on the time variable t for notation convenience, the global, oneparameter semi-discrete DEM scheme takes the form where

Continuous Limit
Due to the substantial disagreement in the scientific community about the governing equations which regulate multiphase phenomena, many authors have tried to derive such mathematical models in different ways. One of the advantages of taking the perspective of the DEM method, is the possibility to derive it, starting from a local description. A first example in this direction was performed in [40], for the specific choice of r ≡ 0. Such a model can be summarized into the following system of PDEs: each phase k = l ∈ {1, 2} evolves according to where the interfacial pressure p I and velocity u I are given by where Z k := ρ k a k denotes the acoustic impedance and mean interfacial pressure p I and velocity u I read Relaxation parameters µ, λ are defined according to the interfacial area In Appendix B we detail the procedure to derive the continuous limit of such scheme, as well as the specific assumptions. The resulting model for the description of (possibly) disperse flow of phase k into l reads Remark 3. An interesting fact concerning this limit is that it has, in the case r = 1, a conservative character, which has already been established by other authors, see [31,38], independently. Indeed, in the limit of small values of α k , one recovers the same model of [31].
In [38], a similar model is proposed replacing the volume fraction equation making assumptions on the production rate of dispersed particles.

Solution Strategy
In this section we make some comments about the resulting scheme (28), its equilibrium variety, its numerical approximation and the use of relaxation procedure. In order to simplify the notation and the following discussion, notice that (28) can be rewritten as i is the numerical contribution coming from the application of the space-discretization operator applied to the states is the average number of internal particles per cell and R is the relaxation term arising from the presence of internal disperse particles. Due to the assumption that the micro-scale is so rich that an infinite number of dispersed particles can be considered inside each cell (i.e. λ i → ∞), the system (34) is typically split into two step, namely the hyperbolic and the relaxation ones. This is also the strategy we follow in this work: the approximation of (28) is accomplished by the following operator splitting method 1. Hyperbolic Step: The hyperbolic step stands for the evolution of the variables according to the left hand side of (34), namely

Relaxation
Step: The relaxation step updates the approximation of the solution U, coming from the hyperbolic step, by computing the equilibrium state of the following ODE as λ i → ∞.
Notice that (35)- (36) are intended also for the volume fraction α k with the formal substitution F ≡ 0 and U ≡ 1.
We conclude this section by stating a convexity property of the numerical scheme (35). In particular, we approximate the set of ODEs (35) with a Forward Euler (FE) method, as it is usual in first-order numerical schemes. Hence, the update formula for (35) reads where we introduced explicitly the dependency of the scheme (28) with respect to the two parameters r n  r). Assume also that each contribution of the numerical flux is positive. Then, the numerical solution (α k U k ) n+1 i predicted by the scheme (35) with the FE time-approximation lies between the corresponding numerical approximations generated by the choices r ≡ 1 and r ≡ 0, i.e.
Proof. Let us rewrite the scheme (28) by using the following form ∈ R is defined in (28). Hence, (37) can be reformulated into By Theorem 3, each P Remark 4. Notice that the above proposition guarantees a bound for the numerical approximation over the hyperbolic step. This is in principle not true for the two-stages scheme (79). Nevertheless, in the following numerical tests we do observe such behavior even though we were not able to prove the conclusion of Proposition 4 when including the relaxation step. This may be due to some monotonicity property of the relaxation step, whose study is out of the scope of this paper.

A comment about the relaxation step
The relaxation step has attracted a lot of attention, due to its paramount importance for an accurate multi-scale description. Due to the discrete nature of the right-hand side in (36), the system of ODEs one has to solve depend on the specific choice of the RS under use. For example, when considering the exact RS (for the single-phase case) one would need to solve two RP (for each computational cell), typically via some root-finding procedure. This latter can become quite cumbersomeness, and a way to circumvent it [28,41,42] is to simplify the system of ODEs (36) by substituting the RS with a fixed, simple approximate RS. Typically the acoustic solver [46,33] constitutes a reasonable and sufficiently simple choice. After such a simplification step, one just derives the corresponding continuous limit of the right hand side (in terms of the variables U i ), and the system is solved by computing the equilibrium state as λ i → ∞. Under the choice of the acoustic solver, one can show that the set of ODEs forces the mixture constituents to move with a single velocity and single pressure, as it was observed/theorized in many works, see [7] and references therein. Unfortunately, there are several simplification steps in this procedure, which do not guarantee that any other reasonable solver leads to the same mechanical effects. This is sometimes reformulated saying that the equilibrium variety (i.e. the set of states U ∞ i at the end of the relaxation step) depends on the choice of the RS. A first investigation in this direction was performed in [4] where the authors showed that for several solvers of common use this is not the case: for such solvers, the equilibrium variety turns out to be defined by the conditions where the index ∞ denotes the states at the end of the relaxation step. Hence, one is tempted to conclude that the relaxation variety is invariant under the choice of (reasonable) solvers.
Here we consider the two following assumptions: 1. Assumption on the Equilibrium Variety: We assume that the equilibrium variety defined by solving (36) and letting λ i → ∞, can be alternatively computed as the reduced set of variables which make R vanish, that is, we assume that there exist a Maxwellian M : 2. Assumption on the Riemann Solver: We assume that the following flux-vector splitting condition holds where Proof. See Appendix C.
Remark 5. Notice that assumption (40) is satisfied by many popular Riemann solvers, including the exact, HLLC, and acoustic solvers.
Remark 6. Notice that (39) does not imply that any solver fulfilling the aforementioned assumptions will produce the same approximations for S ∞ i or p ∞ i . Specifically, different solvers will produce different value for the interface velocities, in general. Thus, the form of the relaxation term is characterizing the equilibrium variety, but it yields no information on how to compute such values.

Numerical Experiments
In this section we test the numerical algorithm to show the influence of the newly derived set of probabilities. Numerical fluxes have been computed using the HLLC flux for the Euler equations [46,47] and Lagrangian fluxes have been computed according to where F HLLC , S * HLLC ,U * HLLC denote the numerical flux, the speed of the contact discontinuity and the intermediate (star) value provided by the HLLC solver, see [46] for details. Notice that the relation (41) is crucial: indeed, one could be tempted to use the acoustic solver provided in [33,40,46] to approximately compute the Lagrangian flux. However, this choice has been found to produce erroneous pressure oscillations near discontinuities, especially in absence of relaxation. Furthermore, we point out that such Riemann Solver for the Lagrangian Flux could also be interpreted to be non-positive conservative in the sense of [16]. For all our simulations we used a CFL constraint of CFL = 0.9. Materials are governed by the stiffened gas equation of state The parameters of gas are γ 1 = 1.4 Pa, π 1 = 0, while for the liquid phase are γ 2 = 4.4 and π l = 6 × 10 8 Pa. Each experiment is computed on the domain D = [−1, 1], unless differently stated.

Uniform Volume Fraction
The first numerical experiment consists of a shock-tube problem, with a uniform volume fraction. The initial mixtures consists of a strong pressure difference. The initial condition in terms of the primitive variables V = [α, ρ, u, p] reads: where ρ 1 = 50, ρ 2 = 1000, p L = 10 9 , p L = 10 5 .

The relaxation-free case
We initially assume that λ i = 0 for any i: solutions associated to a stratified flow evolve independently whereas in the dispersed regime, interactions between the fluids do occur. For the sake of comparison, we report the solution of such problem with r = 0 in Fig.5 (originally proposed in [5]) and the one associated to r = 1 in Fig.6. The latter case is presented using several meshes to show convergence, whereas the case r = 0 is compared to the single-phase exact solutions, to show independence of the two numerical simulations.
The stratified flow regime simulates two non-interacting fluids one on top of the other, while the disperse one models a dilute flow of air inside water. As expected, this latter situation leads to interaction of phases, even though no relaxation is imposed. This is due to the discontinuity of volume fraction at cell interface that enters in the numerical flux through the probability coefficients. Notice the perfect coupling of phases in absence of relaxation for the case r = 1. This clearly highlights the importance of Lagrangian fluxes to maintain it. Furthermore, results for the case r = 1 show near coalescence of velocity and pressure curves: the two phases seem to converge to equilibrium. However, inspection of shock profiles shows slight differences between the fluids, see Fig. 7.
Notice that such an example suggests that, when choosing r = 0, relaxation is not the only mechanical interaction between the two fluids. Finally, an overshoot in the top right corner of gas density phase appears for r = 1, whose amplitude reduces by mesh refinement, see Fig. 7, suggesting convergence in L 1 -norm.

Adding relaxation
We perform the same test, but adding the relaxation procedures described in Appendix D. In this setting an infinite interfacial area is present inside each cell. Results for both regimes (i.e. the stratified and the disperse case) computed with different relaxation procedures are reported in Fig. 8 and    even though discrepancies between the two patterns can be recognized near rarefaction and shocks. Particularly evident is the impact of the choice of r for the post-shock state of density (see first raw of Fig. 9), even if comparable discrepancies can be recognized even for the rest of quantities of interest. Furthermore, discrepancies can also be seen comparing results for different relaxation procedures. For example, differences in shock location predictions are present between relaxation procedures, see Fig.  8. This highlights the fact that the numerical solution of such a test problem is highly dependent on both the relaxation strategy and the probability coefficients at the volume interfaces, raising the question of uniqueness: how can we single-out a physically relevant solution among the infinitely many generated by different realizations of the relaxation strategy and the parameter r ?

Pure Phases
A prototypical benchmark problem for the simulation of two-phase flow is the ability of a scheme of resolving sharp interfaces or reproducing pure phases. Unfortunately the present scheme does not enjoy such property, due to numerical viscosity. Indeed, when attempting to simulate sharp interfaces separating different constituents, the numerical scheme will not maintain the volume fraction in the set X = {0, 1}, due to numerical diffusion. This corresponds to smearing out the interface over several computational cells, thus creating a mixing zone around the exact interface location. A numerical artifact used to circumvent the numerical failure arising in such situation is to assume a negligible amount of dispersed phase, as to stabilize the algorithmic procedure. For the sake of comparison we therefore assume such a strategy to investigate the impact of parameter r when simulating pure phases. We consider the following initial condition in terms of the primitive variables V = [α, ρ, u, p],  mixture quantities of interest are computed for each regime ρ := α 1 ρ 1 + α 2 ρ 2 u := α 1 ρ 1 u 1 + α 2 ρ 2 u 2 ρ p := α 1 p 1 + α 2 p 2 (43) As the sum of the equations for each phase at each point location x = x i results in a formally equivalent system to the single phase Euler equations for the mixture, one can compute the corresponding exact solution according to well-known solvers [22,46]. Results for both the disperse and stratified flow mixtures are compared with exact solution between pure phases in Fig. 10. The continuous-limit relaxation strategy is used for this test case. Again relaxed models produce analogous results, see Fig.  10. Notice that this should be not surprising: indeed, when considering the probability coefficients one these do not depend on the choice of r, in the regions of single-phase flow: let us assume that a material interface is located at x = x i+ 1 2 for some i; then it holds that for any j = i and a k = l ∈ {1, 2} where ε represents the virtual amount of phase l used at the numerical level, and is assumed to be comparatively small, i.e. 0 < ε 1 2 . Then, one gets The significance of this is that, when considering the ideal case of pure phases (ε = 0), one obtains that the only non-zero probability coefficient is the one associated to the probability of having the same phase on both sides of an interface for the phase that has the higher-volume fraction. This corresponds to making all the twophase-fluxes contributions vanish, and the classical, single-phase Godunov scheme is recovered. Hence, each phase behaves independently of the complementary one. Conversely, around the material interface located at x = x i+ 1 2 , it holds α k i = 1 − ε and α k i+1 = ε and α l i = ε and α l i+1 = 1 − ε so that which again make vanish all the contributions not associating to finding phase k and phase l, respectively, on each side of the interface, as one would expect. Notice that such behavior is immediately broken if even negligible (but not-zero) amount of complementary phase is considered in each volume (i.e. ε > 0). This introduces the contribution of other fluxes terms which slightly affect the solution profile, depending on r. Indeed, slight discrepancies can be seen around shocks, even if the overall performance of both models (r = 0 and r = 1) results acceptable and virtually equal. The significance of the above analysis is that, the discrepancies between the two models generated by different choices of parameter r, are dependent on the amount of virtual phase we allocate in each pure chamber. This in turn, also highlights the importance of moderately small disperse particles/sub-scale phenomena in determining shock profiles and corresponding jump relations.

Dynamical creation of interfaces
In this test we examine the capability of the one parameter model to dynamically create interfaces. We consider the following initial condition in terms of the primitive variables V = [α, ρ, u, p], Results for the cavitation test case are reported in Fig. 11, with magnified details shown in Fig. 12.
Here we present results only for the first relaxation procedure of Appendix D. Both the stratified (r = 0) and the disperse (r = 1) regimes are able to dynamically create interfaces, meaning that gas pockets are generated at the discontinuity position. Discrepancies in the velocity field can be appreciated around the discontinuity, like in the oscillations around the peaks of volume fractions. It is worth highlighting that this test presents a moderate speed on both sides of the diaphragm. Increasing the expansion velocity (up to u = 100 m/s, for example) would results in computational failure. Indeed, in the present formulation no mass transfer is considered, so that the creation of gas pockets is only due to the relaxation step.

Randomly chosen, spatially dependent regimes
In this test we want to investigate the difference of predictions with respect to the imposition of randomly chosen r, and a piece-wise constant r. We perform test 1, considering uniform volume  Fig. 13 fraction with first relaxation procedure and compare the results obtained with r = 0 and r = 1 with two constant, randomly chosen r and the following piece-wise constant function Results are shown in Fig. 13: we report only details of the quantity of interest to help appreciate differences. This test yields numerical evidence to understand the impact of the choice of the parameter r on the corresponding numerical approximations. Firstly, one can recognize that solutions display smooth transition as r is increased, when constant throughout space. Notice that the same conclusion carries directly to the piece-wise constant function: depending on the domain of interest, the solution generated by the piecewise constant r lies between the ones computed with constant values, underlying the local dependency of the corresponding solutions.

Dense-to-dilute transition
In this test we aim at investigating the dependency of solutions with respect to variations in the parameter r. Indeed, so far, only spatially constant cases of the parameter r have been considered; here we extend such results to space-time varying functions. Firstly, for the sake of comparison, we fix the same initial condition of Test 1 (i.e. the uniform volume fraction test case), imposing initial stratified flow (r = 0). For each subsequent time level t n , each of the interfacial regime is modeled updating r n for a sufficiently small . Once all the newly generated r n i+ 1 2 are computed, one updates the solution by utilizing the strategy designed in Section 4.2. The parameter encodes the rate of dense-to-dilute transition: negligible values of the parameter will produce virtually same results of the stratified flow case, whereas excessively large values will produce discontinuous flow transition. We present in Fig. 14 only magnified regions of the approximate solutions: results are analogous to the ones reported in Fig. 10, with oscillations at post-shock states. Notice how moderately small values of ε do yield virtually same results as constant r ≡ 0. As suggested by Proposition 4, a convergence towards the stratified flow (r ≡ 0) can be appreciated as decreases, hence demonstrating the smooth dependency at any time of computed solutions with respect to the parameter r. However, the most interesting result of this test is the oscillatory effect appearing for sufficiently large values of . In order to further investigate such oscillatory effect appearing near discontinuities, we run the same test, increasing and comparing results with a uniformly randomly chosen r n i+ 1 2 . Corresponding results are shown in Fig. 15. Oscillatory effects appear near discontinuities which propagate to adjacent regions. As expected, increasing the value of ε will produce with higher probabilities results oscillating around r ≡ 1. Corresponding numerical solutions show indeed variations around the mentioned value r ≡ 1.

Discussion
The modeling of multi-phase flow is very challenging, given the range of scales as well as type of distinct flow regimes that one encounters in this context. We revisit the discrete equation method (DEM) for two-phase flow in the absence of heat conduction and mass transfer. As DEM is based on an ensemble averaging of flow realizations, the mean flow has the potential to describe different two-phase flow regimes. Our starting point was the derivation of Abgrall and Saurel [5] where the authors proposed a DEM for two-phase flow. Our main contributions in this paper was to carefully analyze the resulting probability coefficients and to prove local convexity for them. This rigorously establishes that this version of DEM can indeed model different flow regimes ranging from the disperse to stratified (or separated) flow. Moreover, we reformulated the the resulting mesoscopic model in terms of an one-parameter family of PDEs that interpolates between different regime. The limit cases of this parameter correspond to disperse and stratified flow, respectively. Furthermore, two sets of relaxation procedures were also proposed to enforce relaxation to equilibrium.
We presented extensive numerical experiments to describe the capabilities as well as limitations of the proposed DEM. First, we demonstrated that different values of the probability coefficients yield different predictions on the mechanical interaction between phases. Indeed, it is the probability coefficients, rather than relaxation terms, that lead to this behavior, rather than the details of the relaxation terms which serve to enforce thermodynamic constraints. Indeed, the interaction of phases demonstrated through numerical tests without relaxation procedures, suggest mechanical exchange even if no sub-particles are present inside each volume, in contrast to the interpretation of [5]. This point of view clearly brings out the complimentary roles played by the probability coefficients and relaxation terms in DEM.
The proposed formulation also brings out possible limitations of the DEM approach. In particular, we show that an infinite number of possible models can be constructed, resulting in a ill-posed procedure to construct multiphase simulations. Although several works have investigated the mechanical/thermodynamical consistency of the continuous limit associate to stratified flow, proving it to lead to physically meaningful models. However, even under such an ansatz, the DEM method requires the relaxation operator to be added manually assuming either an infinite drag force or an estimate of the interfacial area in each cell. These latter may become problematic to obtain, if possible, without making any assumption on the flow regime. Indeed, among the desiderata for multiphase flow simulations, the avoidance of user-specification of the flow regime is paramount.
We observe that the DEM scheme represents a finer level of description as compared to the continuum theory approach. In turn, such strategy achieves extensive modelling capabilities, ranging from stratified to disperse flows. However, as demonstrated in this paper, such a mesoscopic approach is not yielding a fully-determined system of constitutive equations as neither equilibrium states nor probability coefficients are uniquely defined. Indeed, this is due to the determination of mean flow variables whereas information about the underlying microstructure is lost. Such inherent under-determination does not rendering the model invalid, but it rather requires closure conditions to be supplied. This is equivalent to saying that the microstructure details lost in the passage to the ensemble averages has to be recovered from somewhere, which in the case of the DEM, is embodied in the probability coefficients and in the relaxation terms. In other words, one needs to adapt the free parameters according to the flow topology, but contrarily to other approaches, it is easy to see where new inputs must be supplied.
Moreover, from a mathematical point of view, the underlying probability coefficients in our DEM scheme can be interpreted as two-point correlation measures. The resulting conclusion is that the DEM approach lacks information about correlation measures, at each space-time location. Therefore, it is out belief that the point of view of measures should be preferred over classical weak forms. Indeed, many recent works dealing with numerical approximations of turbulent flow have shown success of weaker notions than the usual distributional sense. Notice that such an approach is in principle also capable to deal with non-conservative products, typically featured by most well-known two-phase models.
Our results also show that the form of the relaxation variety does not depend on the underlying solver for the hyperbolic step, thus suggesting that a characteristic feature of such phenomena is the determination of the speed at which they reach equilibrium.
Such an insight confirms that the essence of multiphase fluids lies in their microstructure, which has to be considered in order to characterize mean flow variables. For this reasons, forthcoming papers aim at including such information in the modeling of multiphase flow.

A Proof of local convexity
We start the proof of our main result Theorem 3 by reporting the following trivial fact: for any a, b ∈ R it holds max(a − b, 0) + min(a, b) = a Furthermore, we will adopt the notation p = q ∈ {1, 2}.

Proof. (of Theorem 3)
We split the proof into several steps.
2. r does not depend on p: We are going to show that the quotients (50)-(53) are in fact nondepending of p, namely the one induced by the choice p = k coincide with the one induced by p = l, for any k = l ∈ {1, 2}. By the equivalence between (50) and (53), it is enough to show that Subtraction of (22b) from (22a), when p = k and q = l, yields Notice that, and Equations (55), (56) and (57) into the left hand side of (54), leads to the desired equality.

B The Continuous Limit
In this section, we aim at deriving a set of PDEs for the simulation of multiphase flow phenomena. This can be achieved by deriving the continuous limit that the set of discrete ODEs 28 is converging to. As discussed for the relaxation term, the convergence of each single term involved in the system of ODEs is solver-dependent, in principle. One possibility to circumvent such difficulty is to fix a specific form of the RS, which allows for computations. We choose the assumption (40). Furthermore, we make also the following simplification: let us assume that the RS under use computes the contact-discontinuity speed σ and pressure p * as follows where Z k = ρ k c k k = L, R denotes the acoustic impedances computed by the solver and c k is an approximation to the sound speed. Specifically, we always assume that the internal energy can be described in terms of the independent variables ρ k and p k , i.e. e k = e k (ρ k , p k ) denotes the EOS, so that the sound speed is denoted as For the case of the acoustic solver (see. [46], page 299-300), we simply get c k = a k . Conversely, for the case of the HLLC solver (see [46,47]) one has that c k = u k − S k , where S k denotes the fasted signal speed on the k-th side. Notice that each of the aforementioned interfacial solvers can be written into the sum of a symmetric part and anti-symmetric part, namely: so that S(L, R) = S(R, L) and AS(L, R) = −AS(L, R).
We split the analysis into several contributions

B.1 Relaxation Terms
Based on the assumption (40), we get that Simple algebraic manipulations by using (58) and the symmetric-antisymmetric splitting show that Plugging these latter into the relaxation form, one concludes By defining the parameters By assuming that the relative number of interfaces E N int ∆x remains bounded as ∆x → 0, the continuous limit for the relaxation term is derived.

B.2 Conservative Terms
The convergence of conservative fluxes is readily provided: by the finite difference approximation, one gets Under the assumption that kinetic-fluctuations may be disregarded (see 8), one gets that so that the conservative terms of the continuous limits are proven.

B.3 Non-Conservative Terms
Inserting the new set of probabilities the boundary terms can be split into Before detailing each term, we introduce the following convenient notation So that also we rewrite the flux-indicators β (p,q) Under such assumption, the disperse term can be rearranged as due to the discussion of previous subsection. We then focus on the convergence of the second term. Under the hypotheses performed for the previous section, we can write Hence, by analogous splitting to the one performed above, one gets that where α disp denotes the volume fraction with lowest value.

C A solver-invariant equilibrium variety
In this section we are concerned with the proof of a result concerning the equilibrium variety of (36). For the sake of simplicity we will avoid the subscript i, meaning that all the following considerations hold cell-wise. This means U p := (U p ) i for each p ∈ {1, 2}. Furthermore, we will make use of the following notation: F * (U L , U R), σ(U L , U R ), U * (U L , U R ) denote the flux, the interface/contact discontinuity speed and the solution (i.e. the Godunov state) generated from the resolution of the RP by means of a prescribed RS. For the sake of brevity, we will also use q * pq to denote the resulting quantity q in the star region yielded by the resolution of the RP with initial data U L = U p and U R = U q as referring to the RP between the phases p and q. Without loss of generality, we consider the relaxation term (36) for phase k = l ∈ {1, 2}, which reads where the sub-index I denotes the evaluation of the corresponding quantity close to the interface from the side of phase k. By the evaluation of each Lagrangian flux to the interface and assumption (40), one has that Hence, plugging this latter into (63) one gets that the equilibrium variety is defined by the set of ODEs Solving such system of ODEs implies the well-known conditions on relaxed states σ kl = σ lk = S ∞ p * lk = p * kl = p * . Notice that, the first equation in (65) is actually a trivial equation, 0 = 0. This is indeed stating that no matter the values of U ∞ , conditions (63) for mass are always fulfilled. From the point of view of our ODE (36) this implies the following over the relaxation step. Therefore, by assumption on the equilibrium variety, the primitive variables vector of relaxed states can be rewritten as so that the relaxed volume fractions are given by α ∞ k = α 0 k ρ 0 k /ρ ∞ k , by (66). Therefore, a natural Maxwellian M is defined as

D Relaxation Strategies
One important feature of two-phase flow models is to correctly model the interaction between mixture phases. This has been studied for example in [7,36]. Several strategies have been developed so far, and one robust approach is to model interaction between phases by means of relaxation procedure, typically involving stiff source terms. As firstly suggested by Abgrall and Saurel in [5], if the relaxation term E relax F lag i in (28) consists of moderate amount of bubbles, standard resolution of (28) can be applied. However, it is usual to associate such relaxation terms to a large values of source terms, i.e. large numbers of disperse particles are considered. Therefore, relaxation strategies that capture the equilibrium states have to be derived. For the seven-equation model, standard techniques are given in [36,41,28], in which velocity and pressure relaxation steps are split into subsequent operators.
In this work we propose two relaxation strategies that aims at deriving equilibrium states avoiding further splitting methods. 1 ρ k dt. Following the work of [39], a possible choice that has been shown to be compatible with the entropy inequality and with energy conservation is p I (t) ≈ p I (t * ) = p * , u I (t) ≈ u I (t * ) = u * By means of such an approximation, we are led to compute the root of the following non-linear function where e k = e k (ρ k , p) and u * is computed according to (71). The multivariate function F = (F 1 , F 2 ) depends on 2 + 1 = 3 variables, namely ρ k and p, so one equation is missing. We complete the system by enforcing fulfillment of the saturation condition: where, by virtue of mass conservation, α k = (α k ρ k ) 0 ρ k . Hence, the Jacobian matrix of F = (F 1 , F 2 , F 3 ) T reads with definitions The computation of the root of the the multivariate function F is accomplished by means of a standard Newton-Raphson method. The iterative scheme is stopped when the relative increment is sufficiently small and a robust initial guess has been shown to be F 0 = (ρ 10 , ρ 20 , p I (t 0 )) T . Therefore, the approximation of the equilibrium states (ρ * 1 , ρ * 2 , p * ) can be summarized into the following algorithm:

D.2 A projection-relaxation strategy
We propose a second relaxation strategy, that aims at avoiding the computation of the continuous limit of the source term R in (36). This is accomplished making use of the approach developed by Murrone et al in [33]. Such a strategy starts with the introduction of a small parameter to model the speed of the relaxation. More precisely, we introduce the relaxation time ε → 0, so that (34) becomes As discussed in Appendix C, we define a set of relaxed states u which, upon mapping to conserved variables U ∞ = M U (u), defines a root of the function R, namely R M U (u) = 0. Based on the discussion exposed in Appendix C, there exist a natural parametrization in terms of the primitive variables, namely Then, looking for a solution of the form W = M W (u)+ Y, one assumes that there exists an expansion of the source term R such that where J is the Jacobian of the source term R in terms of the variables W evaluated at W = M W (u). Then (79) becomes d dt α If we are able to find the projection matrix P U onto the ker J U R(M (u)), then neglecting the second order terms we get Equations (83) tell us the following: advancing the solution with the hyperbolic step followed by the multiplication of P U is yielding relaxed states. Notice that, this strategy has in principle just the cost of a matrix vector multiplication, in contrast to the rich variety of iterative processes that can arise form solving the continuous limit (36) as λ i → ∞ . One can show that each solver that admits the splitting form (40) are associated to the same projection matrix Π proposed in [4], in case of transonic flow regimes. Indeed, the Jacobian matrix in terms of the primitive variables of M reads We then deduce clearly that the range of J U R(M (u)) is spanned by the vectors However, notice that the Jacobian of the Maxwellian is written in terms of the primitive variables, therefore we need to transform V j in terms of the primitive variables. This can be done, by computing the linear transformation T between conservative and primitive variables (see [4]), such that straightforward computations lead to Assembling the matrix S = [dM (u), T (M (u))V 1 , , T (M (u))V 2 ] and inverting it, yields the projection matrix