A Hierarchy of Non-Equilibrium Two-Phase Flow Models

We consider a hierarchy of relaxation models for two-phase flow. The models are derived from the non-equilibrium Baer-Nunziato model, which is endowed with relaxation source terms to drive it towards equilibrium. The source terms cause transfer of volume, heat, mass and momentum due to differences between the phases in pressure, temperature, chemical potential and velocity, respectively. The subcharacteristic condition is closely linked to the stability of such relaxation systems, and in the context of two-phase flow models, it implies that the sound speed of an equilibrium system can never exceed that of the relaxation system. Here, previous work by Fl{\aa}tten and Lund [Math. Models Methods Appl. Sci., 21 (12), 2011, 2379--2407] and Lund [SIAM J. Appl. Math. 72, 2012, 1713--1741] is extended to encompass two-fluid models, i.e. models with separately governed velocities for the two phases. Each remaining model in the hierarchy is derived, and analytical expressions for the sound speeds are presented. Given only physically fundamental assumptions, the subcharacteristic condition is shown to be satisfied in the entire hierarchy, either in a weak or in a strong sense.


Introduction
The concurrent flow of two fluid phases occurs in a wide range of industrially relevant settings, including in nuclear reactors [1,2], petroleum production [3], heat exchangers [4], cavitating flows [5], and within carbon capture, transport and storage (CCS) [6,7,8]. However, for most simulation purposes, resolving the full three-dimensional flow field may be too cumbersome, due to the complex interaction between the phases. In particular, this encompasses calculating the temporal evolution of the interface between the phases, and the transfer of mass, heat and momentum across it. Averaging methods (see e.g. Ishii and Hibiki [9] or Drew and Passman [10]) may therefore be applied to avoid direct computation of the interface. The resulting coarse-grained models may often be expressed as hyperbolic relaxation systems with source terms accounting for the interactions between the phases, driving them asymptotically towards equilibrium at a finite rate. In a quasi-linear form, one-dimensional versions of such systems may be written as • The non-stiff limit, corresponding to the limit → ∞. In this case, we may write (1) as We will refer to (2) as the homogeneous system.
• The formal equilibrium limit, which is characterized by Q(U) ≡ 0. This defines an equilibrium manifold [13] through M = {U ∈ G : Q(U) = 0}. We now assume that the reduced vector of variables u(x, t) ∈ R n , where n ≤ N, uniquely defines an equilibrium value U = E(u) ∈ M. We may then express (1) as where B(u) = P(u)A(E(u))∂ u E(u) is the Jacobian of the reduced system. Herein, we have defined the operator P(u) : R N → R n through P(u)∂ u E(u) = I n , i.e. the identity matrix. We will refer to (3) as the equilibrium system.
We expect solutions of (1) to approach solutions of (3) as → 0, i.e. in the stiff limit, where the relaxation towards equilibrium is assumed to be instantaneous.

The subcharacteristic condition
An essential concept which arises in the study of relaxation systems and their stability, is the so-called subcharacteristic condition. It was first introduced by Leray [14], subsequently independently found by Whitham [15], and later developed by Liu [16] for conservative 2 × 2 systems. For more general systems, Chen et al. [13] defined an entropy condition which they showed implies the subcharacteristic condition. Yong [17] proved that for n = N − 1, the subcaracteristic condition is necessary for the linear stability of the equilibrium system. Solem et al. [12] proved that it is also sufficient. Hence, for rank 1 relaxation processes, the subcharacteristic condition is equivalent to linear stability.
For a general N × N relaxation system, such as (1), the condition may be formulated as follows.
Definition 1 (Subcharacteristic condition). Let the N eigenvalues of the matrix A of the homogeneous system (2) be given by Λ i , sorted in ascending order as Similarly, let λ j be the n eigenvalues of the matrix B of the equilibrium system (3). Herein, the homogeneous system (2) is applied to a local equilibrium state U = E(u), such that Λ i = Λ i (E(u)), and λ j = λ(u). Now, the equilibrium system (3) is said to satisfy the subcharacteristic condition with respect to the homogeneous system (2) when (i) all λ j are real, and (ii) if the λ j are sorted in ascending order as then the eigenvalues of the equilibrium system are interlaced with the eigenvalues of the homogeneous system, in the sense that λ j ∈ [Λ j , Λ j+N−n ].
The subcharacteristic condition has been shown to be an important trait of many physical models [18,19,20], since the eigenvalues then have a direct physical interpretation as the characteristic wave speeds of the system. In the context of relaxation models for two-phase flow, the fastest wave speeds are the speeds of pressure waves, which are identified as the fluid-mechanical speeds of sound. The subcharacteristic condition then implies in particular that the sound speeds of an equilibrium model can never exceed that of the relaxation model it is derived from. This is consistent with the folklore knowledge in the fluid dynamics community that the "frozen" speed of sound is higher than the equilibrium speed of sound [21,22,23].

The model hierarchy
In a general averaged two-phase flow model, the two-phase mixture will consist of two fluids which evolve independently. We assume local thermodynamic equilibrium in each phase, i.e. that each of the phases may be described by an equation of state, such that specifying two thermodynamic quantities completely determines all thermodynamic properties of that phase. Herein lies also the assumption that the thermodynamic quantities are unaffected by the local velocity field. Each phase k may then be thought of as having separate pressures p k , temperatures T k , chemical potentials µ k , and velocities v k . Since the two-phase mixture will move towards phase equilibrium in each of the mentioned variables, we may model these interactions by employing relaxation source terms corresponding to the following relaxation processes: p volume transfer, i.e. relaxation towards mechanical equilibrium due to pressure differences between the phases, i.e. expansion or compression; T heat transfer, i.e. relaxation towards thermal equilibrium, due to temperature differences between the phases; µ mass transfer: relaxation towards chemical equilibrium due to differences between the phases in chemical potential; 2 and v momentum transfer, i.e. relaxation towards velocity equilibrium, due to velocity differences between the phases, i.e. friction.
The starting point of the forthcoming analysis will be the celebrated Baer-Nunziato (BN) model [24], which is the most general available formulation of a two-fluid model, i.e. models where the phases have separately governed velocities. The BN model is endowed with appropriate relaxation terms corresponding to each of these processes presented above. By considering the homogeneous and equilibrium limits of each relaxation process, i.e. assuming all combinations of zero or more of them to be instantaneous, we obtain a hierarchy of models, each with partial equilibrium in one or more of the aforementioned variables.
This hierarchy can be represented as a four-dimensional hypercube, as illustrated in Figure 1. Herein, each model  Figure 1: The 4-dimensional hypercube representing the model hierarchy. Parallel edges correspond to the same relaxation processes, and each vertex signifies a unique model in the hierarchy, assuming instantaneous relaxation in zero or more of the variables p (pressure), T (temperature), µ (chemical potential) and v (velocity). The leftmost, red circle denoted by "0" represents the Baer-Nunziato model [24]. The colored edges represent relaxation processes where a subcharacteristic condition has been established. This was done by Flåtten and Lund [23] and Lund [25] between the models shown in yellow circles, by Ferrer et al. [26] for the models in green circles, and for the model in blue by Morin and Flåtten [27].
is symbolized by a circle, and corresponds to a "corner" of the hypercube. Parallel edges, in turn, correspond to the same instantaneous relaxation assumption, in the direction of the arrow. The basic model, denoted by "0" and shown in red as the leftmost circle of Figure 1, is thus reducible to all models in the hierarchy. Many of the models in the hierarchy have already been derived, explicitly expressed and thoroughly analyzed, and in this respect, the current paper builds heavily on previous work [28,29,30,31,5,32,33,34]. The models shown in yellow circles in Figure 1 constitute the v-branch of the hierarchy, i.e. the homogeneous flow models, wherein the phase velocities are equal. Such models are subclass of the so-called drift-flux models, where the phasic velocities are related by an algebraic expression. Herein, the v-model was derived by Saurel et al. [34], the vp-model is due to Kapila et al. [30] (see also Refs. [35,29]), and the vpT -model was studied e.g. by [20]. The vpT µ-model is known as the homogeneous equilibrium model and has been studied by several authors, see e.g. Refs. [36,37,38,39,40,41,42]. Flåtten and Lund [23] collected results on the v-, vp-, vpT -, and vpT µ-models, derived the vpµ-model, and showed that the subcharacteristic condition was satisfied for all relaxation processes within this branch of the hierarchy. Lund [25] completed the v-hierarchy by deriving the vT -, vµ-and vT µ-models, and established the subcharacteristic condition in the remainder of the v-branch, given only physically fundamental assumptions.
With regards to the two-fluid models in the hierarchy, several of these models have been derived, employed in simulation [1,3], and analyzed. Here, the p-model was analyzed e.g. in Refs. [43,36], and the pT -model was studied e.g. in Refs. [26,44].
An important issue with p-relaxed (one-pressure) two-fluid models is that they develop complex eigenvalues when the velocity difference between the phases exceeds a critical value, i.e. they become non-hyperbolic [36,45,46,47,27]. This may lead to the lack of stable mathematical and numerical solutions. Nevertheless, these models are extensively used for practical applications; and in numerical simulations they are often mitigated by specifying a regularizing interfacial pressure (see [48,1,49]). Further, estimates of fluid-mechanical sound speeds is of practical importance for the construction of efficient numerical schemes [50,44]. For relations between two-fluid models, we find, as in Ref. [26], the need to state a weaker formulation of the subcharacteristic condition.
Definition 2 (Weak subcharacteristic condition). When the subcharacteristic condition of Definition 1 holds with the additonal equilibrium condition of equal phasic velocities, the weak subcharacteristic condition is said to be satisfied.
The p-and pT -models were analyzed by Martínez Ferrer et al. [26], who showed that the subcharacteristic condition, in a weak or strong sense, is satisfied with respect to existing neighbouring models. Similarly, Morin and Flåtten [27] derived the pT µ-model, and showed that subcharacteristic conditions were satisfied in relation to existing neighbouring models. The highlighted edges in Figure 1 summarize the relations between models where a subcharacteristic condition has already been shown to be satisfied.

Contributions of this paper
The objective of the current paper is complete the study of the subcharacteristic condition in the full hierarchy of two-phase flow models, proving the remaining subcharacteristic conditions. In this respect, a generalization of the work by Flåtten and Lund [23] and Lund [25] is provided, extending the hierarchy to encompass also two-fluid models, i.e. models with separate momentum equations for the two phases. To this end, the derivations of the two-fluid T -, µ-, pµ-and T µ-models represent original contributions. Expressions for the sound speeds in these models are provided. Moreover, we show that the remaining 15 subcharactistic conditions are satisfied, i.e. that the subcharacteristic condition is everywhere respected in the hierarchy, either in a strong or in a weak sense. This is done by comparing the new expressions for the sound speeds to many known results from the literature, and by using techniques involving writing the difference of wave velocities as sums of squares (cf. [25,23]). We present each of the models for which we prove at least one subcharacteristic condition.

Outline
The organization of the current paper is as follows. In Section 2 we present the basic model with all possible source terms, derive evolution equations for the primitive variables, and state a parameter set which suffices to satisfy the laws of thermodynamics. In Sections 3 to 8, we present in turn the v-, p-, T -, µ-, pµ-and T µ-models. For each model we give explicit analytic expressions for the sound speeds, and prove the remaining subcharacterisic conditions with respect to related models. In Section 9 we show plots of the sound speeds in the different models, and briefly discuss physical and mathematical properties of models in the hierarchy. Finally, in Section 10, we draw conclusions and suggest possible future work.

Basic model
In this section, we present the basic BN model [24]. In this model, which is hyperbolic, the two phases have separate pressures, temperatures, chemical potentials and velocities. We state the model in a form reminiscent of that proposed by Saurel and Abgrall [51], but with all four possible relaxation source terms accounting for the interaction between the phases. From this, we determine the evolution equations of the primitive variables. Based on the evolution equations, we derive a parameter set which suffices for the model to satisfy fundamental physical laws.

Governing equations
In the following, we present the governing equations in the basic model, supplemented with physically appropriate relaxation terms. We let α k denote volume fraction, v k velocity, ρ k density, p k pressure, T k temperature, µ k chemical potential, e k internal energy per mass, for each phase k ∈ {g, }, where g denotes gas and denotes liquid.

Volume advection
We assume that apart from advection, the interface between the phases can only move due to pressure differences. This is commonly formulated as wherein v i is an interface velocity and I is the pressure relaxation parameter. Hence, the volume fraction is advected with the velocity v i . There are several discussions available on how to choose this interface velocity, see e.g. [32,52].
In the following, we shall motivate it from a thermodynamic point of view, using the second law of thermodynamics. The local volume transfer must occur so that the phase with the lowest pressure is compressed, and the phase with the highest pressure is expanded. This is enforced through I ≥ 0. Moreover, the volume fractions must satisfy α g + α = 1, where α k ∈ (0, 1), and hence only one evolution equation for the volume fractions is needed.

Mass balance
The evolution of the mass of each phase is contained in the balance equations wherein K is the mass relaxation parameter, and the source terms on the right hand sides of (7) and (8) account for mass transfer between the phases [53,54]. The mass transfer occurs from the phase with the highest chemical potential towards the phase with the lowest, which is ensured through the assumption K ≥ 0. We observe that conservation of total mass is contained by summing (7) and (8): wherein we have defined the mixture density ρ = α g ρ g + α ρ .

Momentum balance
Similar balance laws apply for the momentum of each phase: Herein, p i is an interface pressure and M is the momentum relaxation parameter. Note that from the averaging procedure resulting in these models, the interface velocity v i in (10) and (11) need not be the same as that in (6) (see e.g. Ref. [55]). However, we have chosen these to be equal to keep the notation to a minimum, as this will not influence the main conclusions of this paper. The source terms associated with v i on the right hand sides of (10) and (11) represent the momentum of the condensating or vaporizing fluid, which is transferred to the other phase. The source terms associated with M represent interfacial friction, and are assumed to cause momentum transfer from the phase with highest velocity towards the one with lowest velocity, which is ensured by requiring M ≥ 0. We observe that conservation of total momentum is ensured by summing (10) and (11):

Energy balance
The balance laws for the energy of each phase may be stated as . (14) Herein, µ i is an interface chemical potential, H is the temperature relaxation parameter, and we have introduced the total phasic energy per volume E k = E int k + E kin k , where the phasic internal and kinetic energies are given by, respectively, On the right hand side of (13) and (14), the terms associated with I represent energy transfer due to expansioncompression work, the terms associated with K represent the energy which the condensating or vaporizing fluid brings into the other phase, the terms associated with M represent energy transfer due to frictious momentum transfer, and the terms associated with H represent pure heat flow. The latter should flow from the hotter to the colder phase, which is ensured through the assumption H ≥ 0. Moreover, we see that total energy is conserved by summing (13) and (14), where we have introduced the mixed total energy per volume E = E g + E .

Phase independent form
With all possible relaxation terms, the BN model [24], as presented in (6) to (8), (10), (11), (13) and (14), can be stated compactly as for each phase k ∈ {g, }. Herein, the shorthand forms of the relaxation source terms, I k , K k , H k and M k , have been defined such that

Evolution of primitive variables
In order to systematically derive other models in the hierarchy, and to derive a physically valid parameter set for the basic model, we now seek the evolution equations for primitive variables, such as phasic velocity v k , density ρ k , pressure p k , temperature T k , entropy s k and chemical potential µ k . To simplify the notation in the forthcoming, we introduce the phasic material derivative, defined by for each phase k ∈ {g, }.
In the forthcoming calculations, the following relation will prove useful. For an arbitrary quantity f , we have from (19) and (22) that Note that in order to keep the following analysis as short and concise as possible, most of the details in the derivations are omitted. Further, in the analysis we assume each phase to be described by a completely general (smooth) equation of state. The assumption that any primitive thermodynamic variable can be determined by knowledge of any two other is important (see also the Supplementary Material for a list of relations between thermodynamic differentials). Modelling equations of state and closure laws is beyond the scope of this paper.
We now proceed to find the evolution equations for the primitive variables.

Volume fraction
For clarity we state the evolution equation for the volume fraction. Using (18), we have that

Velocity
We now seek the evolution equation for the phasic velocity. Using f = v k in (23), and (20), we obtain

Density
The density evolution equation is found by combining (19) and (24),

Kinetic energy
In order to obtain the evolution equation for the specific internal energy, we start by finding the evoluton equations for the kinetic energy. Using f = v 2 k /2 in (23), and (16) and (25), we obtain

Internal energy
We obtain the evolution equation for the internal energy by subtracting (27) from (21), expanding and collecting terms: where we have introduced a shorthand interface energy g k = µ i + 1 2 (v i − v k ) 2 . Now, by using (15) and (28) and f = e k in (23), we obtain

Entropy
The fundamental thermodynamic differential-alternatively referred to as the first law of thermodynamics-reads where s k is the specific entropy of phase k. By writing (30) in terms of material derivatives, and inserting (26) and (29), we obtain the evolution equation for the phasic entropy as Herein, the phasic specific enthalpy is defined as h k = e k + p k /ρ k . By using f = s k in (23) along with the identity µ k = h k − T k s k , (31) may be written in the balance form where we have defined the phasic entropy per volume S k = α k ρ k s k .

Pressure
The pressure differential in terms of the density and entropy differentials may be written as where we have introduced the phasic thermodynamic speed of sound and the first Grüneisen coefficient, respectively defined by c 2 By writing (33) in terms of the phasic material derivative, and inserting (26) and (29), we arrive at

Temperature
We now seek the equation governing the phasic temperature evolution. The temperature differential may in terms of the pressure and entropy differentials be written as where the specific isobaric heat capacity is defined by C p,k = T k (∂s k /∂T k ) p k . Now, writing (35) in terms of phasic material derivatives, and inserting (31) and (34), we obtain 2.2.9. Chemical potential The natural differential of the phasic chemical potential reads Therefore, writing (37) in terms of phasic material derivatives, and inserting (34) and (36), we obtain

Laws of thermodynamics
For the model to correctly represent physical phenomena, it should be verified that it satisfies fundamental physical principles [23,20]. We have already verified that it conserves mass, momentum and energy, respectively represented by (9), (12) and (17), where the latter is known as the first law of thermodynamics. We now consider the second law of thermodynamics, which states that the total entropy should be non-decreasing. The analysis in the following is reminiscent of that of various previous works [23,32].

Total entropy evolution
The total entropy per volume is given by S = S g + S . The evolution equation for the total entropy is therefore found by summing (32) over k ∈ {g, }: where we have defined the entropy source terms

The second law of thermodynamics
We define the global entropy as where C ⊆ R is some closed region.
Definition 3. The second law of thermodynamics states that the global entropy is non-decreasing, i.e., in our context. Proposition 1. Sufficient conditions for the relaxation model given by (6) to (8), (10), (11), (13) and (14) to satisfy the second law of thermodynamics (Definition 3) are given only the physically fundamental assumption T k ≥ 0 for k ∈ {g, }.
Proof. By temporal differentiation of (44), in combination with (39) and (45), we obtain where we have assumed that the entropy flux of (39), S g v g + S v , vanishes at the boundary of C . For (50) to be satisfied, clearly S ≥ 0 is a sufficient criterion, for which statement to hold the non-negativity of all the partial source terms S p , S µ , S T and S v is in turn sufficient. We now show this for each of the terms under the conditions of (46) to (49). Firstly, the conditions (48) and (49) inserted into (40) yields Now, (47) is equivalent to µ i = β µ µ g + (1 − β µ )µ , with β µ ∈ [0, 1]. Hence, combination of (41) and (49) yields Next, (49) inserted into (42) yields Finally, (43) becomes and hence all the source terms are non-negative. (48) and (49) are sufficient, not necessary, and the square-root-of-temperature weighted average between the phasic values differs from choices in the literature, e.g. the initial choices by Baer and Nunziato [24]. The reason for this particular weighting is that we enforced the interface velicities in (6), (10) and (11) to be equal, as noted above. Allowing these to differ would enable other linear combinations of the phasic quantities, which could possibly be more suitable for numerical simulations [52]. These differences, however, do not have implications for the main conclusions of this paper.

Wave velocities
We now consider the homogeneous limit of the BN model, where the source terms I , K , M , H → 0. The resulting model has previously been extensively studied by several authors, see e.g. [51,33]. The model has two fluid-mechanical sound speeds; one for each of the phases. The seven wave velocities are given by In typical applications, the flow is subsonic, i.e. |v g − v | c g , c may be a valid approximation. Evaluated in the velocity equilibrium limit, taking v ≡ v g = v , the eigenvalues are, sorted in ascending order, where we have defined c 0,+ = max{c g , c } and c 0,− = min{c g , c } as the higher and lower sound speeds, respectively.

The v-model
We now study the model that arises upon imposing instantaneous equilibrium in velocity, i.e. letting the velocity relaxation parameter M → ∞, which we expect corresponds to Simultaneously, we require the term M g = M (v − v g ) to remain finite. By noting that for a general function f , the phasic material derivatives are equal for the two phases, i.e. D k f = ∂ t f + v∂ x f ≡ D f , then the system that results from evaluating (25) for the two phases k ∈ {g, } can be solved to yield where we have introduced the phasic mass fractions Y k = α k ρ k /ρ. The model that now results from inserting (56) and (57) into the basic model of Section 2, was analyzed by Flåtten and Lund [23,25], as it constitutes the basic model of the v-branch of the hierarchy. The model is hyperbolic and has previously been studied by many authors [56,30,34].

Wave velocities
The wave velocities of the velocity equilibrium model, in the homogeneous limit where I , K , H → 0, are given by [23] λ Herein, the sound speed of this model is defined by Proposition 2. The v-model satisfies the subcharacteristic condition with respect to the basic model, given only the physically fundamental conditions ρ k , c 2 k > 0, for k ∈ {g, }. Proof. We observe that Y g + Y = 1, and due to the given positivity conditions, we have that Y k ∈ (0, 1). Therefore, (59) implies that min{c g , c } ≤ c v ≤ max{c g , c }. It then follows trivially that the wave velocities of the v-model are interlaced in the wave velocities (55) of the basic model evaluated in the velocity equilibrium state (56). Hence, the associated subcharacteristic condition of Definition 1 is satisfied.

The p-model
In this section, we consider the mechanical equilibrium model, which arises when we assume instantaneous mechanical equilibrium in the basic model of Section 2. We let the pressure relaxation parameter I → ∞, which we expect to correspond to p g = p ≡ p. Simultaneously, the product I g = I (p g − p ) should remain finite. The mechanical equilibrium model is found by using (34) evaluated for each of the two phases. From this, we may find an expression for I g without temporal derivatives, and insert it into the basic model of Section 2. The resulting model has been extensively studied previously [51,26]. As other one-pressure two-fluid models, the model is not hyperbolic.

Wave velocities
We consider now the homogeneous limit, where K , M , H → 0. The eigenvalues to the lowest order in the small parameter ε = v g − v , i.e. evaluated in the equilibrium state defined by (56), are given by [26] where the sound speed in the p-model is given by Proposition 3. The p-model satisfies the weak subcharacteristic condition of Definition 2 with respect to the basic model of Section 2, subject only to the physically fundamental conditions ρ k , c 2 k > 0, for k ∈ {g, }, in the equilibrium state defined by (56).
Proof. We see from (61) that c 2 p is a convex combination since ϕ g + ϕ = 1, and ϕ k ∈ (0, 1), due to the given conditions. This implies that and hence the weak subcharacteristic condition is fullfilled with respect to the basic model, whose local eigenvalues evaluated in the same state are given by (55).

The T-model
In this section, we investigate the thermal-equilibrium model (T -model), which emerges from assuming instantaneous thermal equilibrium in the basic model of Section 2. To this end, we let H → ∞ herein, which we expect corresponds to in such a way that H g = H (T − T g ) remains finite. In the following we present the governing equations.

Governing equations
The full T -model may be stated as the basic model of Section 2, in which (13) and (14) are replaced by (17) and the thermal equilibrium condition (64).
In order to establish the impact of instantaneous thermal relaxation on the wave velocities, we need to express the model on a quasi-linear form, and thus obtain the velocities as the eigenvalues of the associated Jacobian. This is most easily done by exploiting the primitive variables, which is what we now turn to do.
Firstly, we have that the phasic pressure differential in terms of density and temperature may be written as where we have introduced the ratio of specific heats ζ k = 1 + Γ 2 k C p,k T/c 2 k , and used (64). With (65), (25) becomes where we have defined the phasic mass per volume m k = α k ρ k , the phasic interface pressure jump ∆ i p k = p i − p k , and the phasic interface velocity difference where we have introduced the extensive heat capacity at constant pressureC p,k = m k C p,k . We now define the weighting factor θ k =C p,k ζ −1 k /(C p,g ζ −1 g +C p, ζ −1 ), for which clearly θ g + θ = 1 and θ k ∈ (0, 1). Multiplying (67) by θ k , and summing over the phases yields where have used the interface parameter definitions of (48) and (49) evaluated in thermal equilibrium (64) to simplify.

Wave velocities
We now seek the wave velocities, i.e. eigenvalues, in the homogeneous limit, where the relaxation source terms I , K , M → 0. From (24), it is then clear that α g is a characteristic variable of the system, since the volume fraction is advected with the velocity v i in the absence of relaxation source terms. By using (26), (66) and (68), the remaining, reduced system may now be expressed on the quasi-linear form ∂ tũT +Ã T (ũ T ) ∂ xũT = 0, whereũ T = [ρ g , ρ , v g , v , T ], and the associated Jacobian is given bỹ from which we can find the remaining five eigenvalues. The characteristic polynomial of the latter is a fifth-degree polynomial, for which in general no closed-form solution can be obtained. We now note that we may writeÃ T = The matrices are given bỹ Hence, we approximate the eigenvalues by means of a perturbation expansion in the small parameter ε. To the lowest order in ε, v g = v =v = v, and the eigenvalues of the T -model are given by where the two distinct sound speeds of the model are given by Proposition 4. The T -model satisfies the weak subcharacteristic condition with respect to the basic model of Section 2, subject only to the physically fundamental conditions ρ k , C p,k , T > 0, for k ∈ {g, }, in the equilibrium state defined by (56).
Proof. We first show that the sound speeds are real. We note that on the given conditions, clearly c 2 T,± ∈ R, and moreover, c 2 T,+ ≥ 0. The product of the sound speeds may be written as Based on the given conditions, it is clear that Z 0 T ≥ 0 and therefore and hence also c 2 T,− ≥ 0, and thus c T,± are real, and by definition, positive. Now, using the definitions of c 0,± and (72), it follows that where The given conditions ensure that Q 0 T ≥ 0. The only ordering of sound speeds compatible with (74) and (75) is 0 ≤ c T,− ≤ c 0,− ≤ c T,+ ≤ c 0,+ , and hence the subcharacteristic condition of Definition 1 is satisfied.
Proposition 5. The vT -model of Lund [25] satisfies the subcharacteristic condition with respect to the T -model, given the physically fundamental assumptions ρ k , C p,k , T > 0, for k ∈ {g, }.
Proof. The sound speed of the vT -model is given by [25] Now, using (72), we can write the product of the differences as where With the given conditions, clearly Q T vT ≥ 0. Hence exactly one of the factors on the left hand side of (78) is negative, and combined with Proposition 4 we realize that c T,− ≤ c vT ≤ c T,+ , and hence the subcharacteristic condition is satisfied.
Proposition 6. The pT -model satisfies the weak subcharacteristic condition with respect to the T -model, given the physically fundamental assumptions ρ k , C p,k , T > 0 for k ∈ {g, } in the equilibrium state defined by (56).
Proof. The sound speed of the pT -model is given by [26] We may now write where Clearly Q T pT ≥ 0, on the given conditions. Hence exactly one factor on the left hand side of (81) is negative, yielding c T,− ≤ c pT ≤ c T,+ , and the weak subcharacteristic condition is satisfied.

The µ-model
We now proceed to investigate the chemical-equilibrium model (the µ-model), which arises when we assume instantaneous chemical equilibrium, i.e. let the chemical relaxation parameter K → ∞, which we expect corresponds to Simultaneously, we require the product K g = K (µ − µ g ) to remain finite, and in the forthcoming we seek to express this without any temporal derivatives.

Remark 2.
It should be noted that there does not seem to be a general agreement in the literature on how to properly model mass transfer (see e.g. [57, pp. 13]). Strictly enforcing (83) may sometimes lead to unphysical results [58]. The present choice (83) is primarily motivated by compliance with the subhierarchy compiled by Lund [25], and evaluating the physical relevance of these models is out of the scope of the present work.
The chemical potential evolution equation (38) may be written as where we have used (83), and defined the shorthands By using (84) evaluated for each of the phases, and subtracting these expressions from each other, we obtain where we have defined the shorthands 6.1. Governing equations By using the expression (86) to insert for K g in the basic model of Section 2, the µ-model can now be summarized with the following set of equations: • Volume advection: ∂ t α g + v i ∂ x α g = I g , • Conservation of mass: ∂ t ρ + ∂ x α g ρ g v g + α ρ v = 0, • Momentum balance: • Energy balance: Momentum and energy equations for the liquid phase are found by phase symmetry; interchanging indices g and .

Evolution of primitive variables
In order to write the system in a quasi-linear form, and thereby find the wave speeds of the µ-model, we use the evolution equations for the primitive variables. We therefore now seek the evolution of some of the primitive variables under the assumption of instantaneous chemical equilibrium.
We first define the weighting factor φ k = χ −1 k /(χ −1 g + χ −1 ). Multiplying (84) by φ k and summing over the phases, we get for the chemical potential For the phasic velocity v g , we find from (25) the evolution equation and v is found by phase symmetry.
The phasic pressure evolution is found from (34). For the gas phase, it reads The corresponding expressions related to the evolution of p are found by phase symmetry.

Wave velocities
We now wish to derive the wave velocities of the µ-model in the homogeneous limit, where I , H , M → 0. In this limit, the volume fraction α g is a characteristic variable with the associated eigenvalue v i . The remaining, reduced model, i.e. (92), (94) and (95) for both phases, may then be expressed in the quasi-linear form ∂ tũµ +Ã µ (ũ µ )∂ xũµ = 0, where the reduced vector of unknowns isũ µ = [µ, v g , v , p g , p ], and the reduced Jacobian reads Again the eigenvalues λ are given the roots of a fifth degree polynomial, for which in general no closed-form solution exists. We therefore expand in the small parameter ε = v g − v , i.e.Ã µ =Ã (0) µ + εÃ (1) µ + . . ., and λ = λ (0) + ελ (1) + . . .. Herein, the lowest-order system matrix reads, takingv = φ g v g + φ v , where we have used the lowest-order term of κ µ , as defined in (87): To the lowest order in ε, v g = v =v = v, and thus the eigenvalue problem consists in finding the roots of det(Ã (0) µ − λ (0) I) = 0. Hence, the full vector of eigenvalues is given by where the two sound speeds in the µ-model are given by Proposition 7. The µ-model satisfies the weak subcharacteristic condition with respect to the basic model of Section 2, given only the physically fundamental conditions ρ k , C p,k , T k > 0 for k ∈ {g, }, in the equilibrium state defined by (56).
Proof. We first note that c 2 µ,± ∈ R on the given conditions, and that c 2 µ,+ ≥ 0. The product of the sound speeds may be written as Given the conditions we have that Z 0 µ ≥ 0, and hence 0 ≤ c 2 µ,+ c 2 µ,− ≤ c 2 0,+ c 2 0,− .
Therefore also c 2 0,− is positive, and thus we have that c 0,± are real and, by choice, positive. Now, the product of the differences of the sound speeds may be written as where Clearly, with the given conditions, Q 0 µ ≥ 0, and hence the only ordering of sound speeds compatible with (106) and (107) is 0 ≤ c µ,− ≤ c 0,− ≤ c µ,+ ≤ c 0,+ , which means that the weak subcharacteristic condition is satisfied.
Proposition 8. The vµ-model satisfies the subcharacteristic condition with respect to the µ-model, subject only to the physically fundamental conditions ρ k , C p,k , T k > 0, for k ∈ {g, }.
Proof. The sound speed in the vµ-model is given by [25] We now consider the product of the differences in the sound speeds of the two models, which may be written as where Clearly Q µ vµ ≥ 0. Hence, exactly one of the factors on the left hand side of (110) must be negative, which gives c µ,− ≤ c vµ ≤ c µ,+ , i.e. the subcharacteristic condition is satisfied.

The pµ-model
We now consider the model which arises when we impose instantaneous mechanical-chemical equilibrium, i.e. we let the relaxation parameters I , K → ∞, which we expect corresponds to p g = p ≡ p and µ g = µ ≡ µ. (112) Simultaneously, I g = I (p g − p ) and K g = K (µ − µ g ) should remain finite. We now seek explicit expressions for these terms in order to find the governing equations of the model. In the following analysis we use the parameter set stated in Section 2 and therefore let the interfacial pressure jump ∆ i p = p i − p = 0. From (34) and (84) we have where we have definedĨ Eqs. (113) and (114) evaluated for each phase now constitute a 4 × 4 system which is straightforward to solve for the four unknowns ∂p/∂t, ∂µ/∂t,Ĩ g , and K g , in terms of spatial derivatives and the remaining source terms. The final expressions for the latter two arẽ where the coefficients are given in Appendix A.

Governing equations
Inserting the expressions (115) and (116) into the basic model of Section 2, we are now in a position to state the full model. The mechanical-chemical equilibrium model may thus be formulated as follows.
• Conservation of mass: • Momentum balance: • Energy balance: The momentum and energy equations for the liquid phase are found by phase symmetry.

Wave velocities
We now wish to write the system in a quasilinear form, in order to find the wave speeds of the system, in the homogeneous limit where we let the relaxation terms M , H → 0. To this end, we will express the model in the vector of unknowns u pµ = [p, µ,v, v g , v ]. We therefore seek the evolution equations for the elements of u pµ .
For the volume evolution, we find, using (24) and (115), that For the volume-averaged velocityv we find, using (25), (115), (116) and (119), that where we have defined the shorthand coefficients P pμ v , G pμ v , V pμ v,g (for which expressions are given in Appendix A), used ε = v g − v , and inserted β g = 1 − β = T /( T g + T ). Now, for the pressure and chemical potentials, we get from (113) and (114) that Again, the coefficients are given in Appendix A.
The homogeneous system in a quasilinear form thus reads ∂ t u pµ + A pµ u pµ ∂ x u pµ = 0, where the system Jacobian is given by Obtaining the assocated eigenvalues exactly by analytic means is again unfeasible, as the problem consists in finding the roots of a fifth-degree polynomial. We therefore expand in ε: A pµ = A (0) pµ + εA (1) pµ + ε 2 A (2) pµ + . . ., where it is assumed that the matrices A (i) pµ are independent of ε. To the lowest order, where ε → 0, taking v = v g = v , we get the matrix where the superscript "(0)" on the coefficients signifies the zeroth-order expansion in ε, such that The eigenvalues in the pµ-model are, to the lowest order in ε, where we have identified the sound speed c pµ of the model, given by Proposition 9. The pµ-model satisfies the weak subcharacteristic condition with respect to the p-model, given only the physically fundamental conditions ρ k , C p,k , T k > 0, for k ∈ {g, }, in the equilibrium state defined by (56).
Proof. From (61) and (127), we observe that we may write Due to the given physical conditions, Z p pµ ≥ 0, and hence 0 ≤ c pµ ≤ c p , i.e. the subcharacteristic condition is satisfied.
Proposition 10. The pµ-model satisfies the weak subcharacteristic condition with respect to the µ-model, under the physically fundamental conditions ρ k , C p,k , T k > 0, for k ∈ {g, }, in the equilibrium state defined by (56).
Proof. Using the expressions (104) and (127) for the sound speeds in the two models, we may write where Clearly, on the given conditions, Q µ pµ ≥ 0. Therefore, exactly one factor on the left hand side of (129) is negative, and hence c µ,− ≤ c pµ ≤ c µ,+ , so the subcharacteristic condition is satisfied.
Proposition 11. The vpµ-model satisfies the subcharacteristic condition with respect to the pµ-model, given the physically fundamental conditions ρ k , C p,k , T k > 0.
Proof. The sound speed in the vpµ-model is given by [23,25] Now, we may write which is clearly positive, due to the given conditions. Thus, 0 ≤ c vpµ ≤ c pµ , i.e. the subcharacteristic condition is satisfied.
Remark 3. By direct comparison of (127) and (131), we find the ratio This is exactly the same ratio as has been shown to hold for other models associated with v-relaxation in the p-branch of the hierarchy [27,26]. We can thus extend the relation by the newly obtained ratio (133) between the sound speeds of the vpµ-and pµ-models.

Proposition 12.
The pT µ-model satisfies the weak subcharacteristic condition with respect to the pµ-model, given the physically fundamental conditions ρ k , C p,k , T > 0, in the equilibrium state defined by (56).
Proof. In the equilibrium state defined by the pT µ-model, we have T g = T ≡ T . The sound velocity in the pT µmodel is given in [27], and may be rewritten as where we have introduced the enthalpy difference ∆h = h g − h .
We may reorganize the last equality in (134) to yield c pµ c pT µ = c vpµ c vpT µ .
Flåtten and Lund [23] showed that the subcharacteristic condition is satisfied between the models on the right hand side, i.e. that 0 ≤ c vpT µ ≤ c vpµ . The same must hold for the models on the left hand side of (136), i.e. 0 ≤ c pT µ ≤ c pµ , and hence the subcharacteristic condition is satisfied. In particular, we may write the sound speed as where Z pµ pT µ =C p,gC p, T 1 ∆h Clearly, Z pµ pT µ ≥ 0 based on the given conditions.

The Tµ-model
We now investigate the model which arises when we assume instantaneous thermal-chemical equilibrium, i.e. let the relaxation parameters K , H → ∞, which expectedly corresponds to T g = T ≡ T and µ g = µ ≡ µ. (139) The products H g = H (T − T g ) and K g = K (µ − µ g ) remain finite, and may be expressed in terms of spatial derivatives and remaining source terms. In the forthcoming, we seek explicit expressions for these terms to insert into the basic model of Section 2.
The equilibrium conditions are contained in (67) and (84). These may be combined to yield where the coefficients are given in Appendix B.

Governing equations
We are now in a position to state the T µ-model in its entirety, by inserting (140) into the basic model of Section 2. The model can be expressed by the following equation set: • Conservation of mass: ∂ t ρ + ∂ x α g ρ g v g + α ρ v = 0, • Conservation of momentum: • Conservation of energy: ∂ t E + ∂ x E g v g + E v + α g v g p g + α v p = 0.

Wave velocities
We now seek the wave velocities of the model in the homogeneous limit, where I , M → 0. As usual, we are interested in the zeroth-order expansion in ε = v g − v . 3 We may therefore directly evaluate the evolution equations in this limit, and take v g = v = v if they are outside the differential operator.
After some tedious, but fairly straightforward algebra, we find that to the lowest order in ε, the wave velocities of the T µ-model are given by Herein, c T µ,± are the sound speeds of this model, which may be expressed by Proposition 13. The T µ-model satisfies the weak subcharacteristic condition with respect to the T -model, given the physically fundamental conditions ρ k , C p,k , T > 0, in the equilibrium state defined by (56).
Proof. We may write where Moreover, we may write where Clearly, Z T T µ ≥ 0 and Q T T µ ≥ 0 on grounds of the given conditions. The only possible ordering of sound speeds is thus 0 ≤ c T µ,− ≤ c T,− ≤ c T µ,+ ≤ c T,+ , i.e. the weak subcharacteristic condition is satisfied.

Proposition 14.
The T µ-model satisfies the weak subcharacteristic condition with respect to the µ-model, given the physically fundamental conditions ρ k , C p,k , T > 0, in the equilibrium state defined by (56).
Proof. We note that we may write where Now, we may also write c 2 µ, Clearly, Q µ T µ ≥ 0 and Z µ T µ ≥ 0 for the given conditions. The only possible ordering of the sound speeds is therefore 0 ≤ c T µ,− ≤ c µ,− ≤ c T µ,+ ≤ c µ,+ , i.e. the eigenvalues of the relaxed model are interlaced between the eigenvalues of the parent model, and the weak subcharacteristic condition is satisfied.

Proposition 15.
The vT µ-model satisfies the subcharacteristic condition with respect to the T µ-model, given the physically fundamental conditions ρ k , C p,k , T > 0.
Proof. The sound speed in the vT µ-model is given by [25] We may now write the product of the differences between the sound speeds as where Due to the given conditions, it is clear that Q T µ vT µ ≥ 0, and thus exactly one of the factors on the left hand side of (152) is negative. Hence, the sound speeds must be ordered as c T µ,− ≤ c vT µ ≤ c T µ,+ , i.e. the subcharacteristic condition is satisfied.
Proposition 16. The pT µ-model satisfies the weak subcharacteristic condition with respect to the T µ-model, subject to the physically fundamental conditions ρ k , C p,k , T > 0, in the equilibrium state defined by (56).
Proof. The sound speed in the pT µ-model is stated in (135). We note that we may write where Q T µ pT µ =C p,gC p, T 2 Due to the given conditions, it is clear that Q T µ pT µ ≥ 0. Hence, exactly one of the factors on the left hand side of (154) must be negative, and therefore c T µ,− ≤ c pT µ ≤ c T µ,+ , i.e. the weak subcharacteristic condition is satisfied.

Comparison and discussion of models
In this section, we compare the models in the hierarchy. We first show plots for relevant cases, and then briefly discuss physical and numerical aspects of the different models.

Comparison of sound speeds
We now present plots of the sound speeds in all the models in the hierarchy, for different physically relevant conditions, in order to illustrate the effect of different equilibrium assumptions. Plots with the same physical parameters were presented by Lund [25] for the v-branch of the hierarchy, building on plots by Flåtten and Lund [23]. Ferrer et al. [26] and Morin and Flåtten [27] presented similar plots for the p-branch of the hierarchy. In the following, results for the complete hierarchy are shown.
The main panel of Figure 2 shows the fluid-mechanical sound speeds of all the models in the hierarchy for a watersteam mixture at atmospheric conditions. The thermophysical parameters are shown in Table 1. Apart from the fact that the subcharacteristic condition is always respected, we notice that there are mainly two equilibrium assumptions that affect the propagation speeds, namely those of p-and v-relaxation, respectively. First, imposing instantaneous equilibrium in pressure attracts the sound velocities to the lower bound of the parent models, which can be seen in the insets of Figure 2. Further imposing velocity equilibrium, the sound velocity is reduced by a factor (see remark 3) The approximation made is valid when ρ g ρ , which is the case here. Hence, these equilibrium assumptions seem to have severe impact on the wave velocities, in particular when the density difference between the phases is large. In Figure 3, the sound speeds in the entire hierarchy are plotted for a CO 2 mixture at 50 bar, whose thermophysical properties are given in Table 2. By close inspection, it can be seen that the subcharacteristic condition is everywhere respected. In particular, the sound speeds of an equilibrium system are always interlaced between the sound speeds in the parent models. Again, the pressure relaxation has the most prominent effect on the sound speed, but also combining thermal and chemical equilibrium seems to have a strong effect.

Discontinuous sound speeds
All the models considered in the present paper are only strictly valid when the gas fraction α g ∈ (0, 1). One would expect the sound speeds of the models to be continuous at the phase boundary, i.e. at the transition between single and two-phase flow, in the sense that the two-phase speed of sound should reduce to the single-phase speed of sound in the limit where one phase disappears: lim for a given model X in the hierarchy. However, some of the models have wave speeds that are discontinuous at the phase boundary. In particular, this concerns the pT µand vpT µ-models, whose sound speeds are discontinuous in both limits α k → 1, which can be seen directly by evaluating the analytic expressions in these limits (see Refs. [25,27]).
The T -and µ-models have "half-continuous" sound speeds, in the sense that for the "±" sound waves, only one of them is continuous in the limit α k → 1. For the µ-model, taking α → 1 in (104) yields lim α →1 c 2 µ,± = T g s 2 g C p,g (c 2 g + c 2 ) + ξ 2 g c 2 c 2 g ± T g s 2 g C p,g (c 2 g − c 2 ) − ξ 2 g c 2 c 2 g 2 T g s 2 which is equivalent to lim α →1 , and lim α →1 Clearly, only one of these approach the appropriate phasic value c 2 . The result limits for α g → 1 are found by phase symmetry. Similarly, we find for the T -model, from (72), that lim α →1 c T,+ = min to which the same observation applies.
The remaining sound speeds are continuous at the phase boundary; for the T µ-model in the sense that lim α k →1 c T µ,+ = c k and lim α k →1 c T µ,− = 0, which can be deduced from the analytic expression (143).

Physical considerations
It is commonly argued that pressure relaxation is a much faster process than the other relaxation processes [5,33]. Temperature relaxation, or heat flow, is associated with diffusion, which is an intrinsically slow process. Chemical potential relaxation, i.e. mass transfer, is also slow compared to pressure relaxation. Zein et al. [33] provide interesting discussions on the topic and argue that temperature relaxation is faster than chemical relaxation. Generally, the magnitudes of the different relaxation times may be strongly problem-dependent. Such considerations may have implications, e.g., for the mass flow through a nozzle, which has been shown to be linked to the subcharacteristic condition [59].
Apart from this, effects not captured by the coarse-grained flow models may come into play, and which model is more accurate may depend heavily on the flow regime under consideration. The effects that arise from having independent phasic pressures may be of importance for the wave dynamics of the system, and thus models with different pressures may be sensible, even though the associated relaxation time is commonly thought to be comparatively short. With regards to evaluating the physical relevance of the models presented herein, experimental data on sound speeds in two-phase flow can be found for various systems [60,61,62,63].

Numerical considerations
A well-known problem with p-relaxed (one-pressure) two-fluid models is that they develop complex eigenvalues when v g v . This is commonly resolved e.g. by adding a regularising pressure which enforces hyperbolicity [1,43,64,65,49]. It is worth noting that the two-fluid models with independent phasic pressures, i.e. the T -, µand T µ-models, are locally hyperbolic even for small perturbations away from velocity equilibrium, due to the following argument: An eigenvalue of a matrix with real coefficients may only be complex if its complex conjugate is also an eigenvalue. Since the eigenvalues of the individual phasic pressure models are real and distinct when ε = v g − v = 0, they must remain so for sufficiently small ε, as the eigenvalues may only become complex in a continuous way. In order to determine how large ε may be before hyperbolicity is lost, we must find the higher-order corrections in ε to the eigenvalues, which is beyond the scope of this work.

Conclusions and further work
In this paper, we have presented and completed a hierarchy of relaxation models for two-phase flow, which arises when we impose instantaneous equilibrium in different combinations of velocity, pressure, temperature and chemical potential. The starting point of the analysis has been the classic seven-equation Baer-Nunziato model [24] equipped with relaxation source terms. We have in the present work provided the T -, µ-, pµ-, and T µ-models, which represent original contributions to the hierarchy. Explicit expressions for the sound speeds of these models have been derived. Using the new expressions and results from the literature, we have shown analytically that the subcharacteristic condition is always satisfied in the hierarchy, given velocity equilibrium between the phases. To this end, we have contributed with 15 new subcharacteristic conditions, stated in propositions 2-16. Out of these, five have been shown in a strong sense, and nine hold in a weak sense, i.e. given equilibrium in velocity.
In further work, the hierarchy could be extended to multi-component or multi-phase flow. Moreover, the different models could be implemented and studied numerically for relevant cases (cf. [52]). Upon comparison with experimental data, one may then unravel under which conditions the different models are admissible. I pµ v = α g α ρ g c 2 g ρ c 2 κ pµ s g C p,g ξ g α g + ξ α + Γ g α g