A hyperbolic two-fluid model for compressible flows with arbitrary material-density ratios

Abstract A hyperbolic two-fluid model for gas–particle flow derived using the Boltzmann–Enskog kinetic theory is generalized to include added mass. In place of the virtual-mass force, to guarantee indifference to an accelerating frame of reference, the added mass is included in the mass, momentum and energy balances for the particle phase, augmented to include the portion of the particle wake moving with the particle velocity. The resulting compressible two-fluid model contains seven balance equations (mass, momentum and energy for each phase, plus added mass) and employs a stiffened-gas model for the equation of state for the fluid. Using Sturm's theorem, the model is shown to be globally hyperbolic for arbitrary ratios of the material densities $Z = \rho _f / \rho _p$ (where $\rho _f$ and $\rho _p$ are the fluid and particle material densities, respectively). An eight-equation extension to include the pseudo-turbulent kinetic energy (PTKE) in the fluid phase is also proposed; however, PTKE has no effect on hyperbolicity. In addition to the added mass, the key physics needed to ensure hyperbolicity for arbitrary $Z$ is a fluid-mediated contribution to the particle-phase pressure tensor that is taken to be proportional to the volume fraction of the added mass. A numerical solver for hyperbolic equations is developed for the one-dimensional model, and numerical examples are employed to illustrate the behaviour of solutions to Riemann problems for different material-density ratios. The relation between the proposed two-fluid model and prior work on effective-field models is discussed, as well as possible extensions to include viscous stresses and the formulation of the model in the limit of an incompressible continuous phase.


Introduction
The difficulties in developing hyperbolic two-fluid models for disperse multiphase flows has been reviewed by Lhuillier, Chang & Theofanous (2013). Many of the models that have been proposed in the literature suffer from being mathematically ill posed (see Drew & Passman 1998;Vazquez-Gonzalez, Llor & Fochesato 2016, for other discussions of this topic), most notably when the Archimedes force is included. Mathematically, well posedness of nonlinear multiphase flow models implies hyperbolicity of the underlying Cauchy problem (Métivier 2005). In practice, numerical simulations with non-hyperbolic two-fluid models diverge under grid refinement due to the complex eigenvalues in the continuum limit (see, e.g. Ndjinga 2007;Kumbaro & Ndjinga 2011). To solve this problem, ad hoc correction terms have been added to make the models well posed (see, e.g. Panicker, Passalacqua & Fox 2018). In particular, some authors have resorted to neglecting the Archimedes force (see, e.g. Hank, Saurel & Le Metayer 2011), which is the root cause of non-hyperbolicity. For bubbly flows the Archimedes force is of critical importance when buoyancy effects are present.
Starting from a kinetic-theory description, Fox (2019) developed a hyperbolic two-fluid model for gas-particle flows that neglects added-mass effects (as well as inelastic collisions and viscous effects). The model equations were derived starting from the Boltzmann-Enskog kinetic theory for a binary hard-sphere mixture. A closure for the particle-pair distribution functions was introduced to account for the Archimedes force in the limit where one particle diameter is much smaller than the other. However, because the closure for the particle-pair distribution function only accounts for mean gradients, it cannot capture the higher-order correlations needed for added mass. The system of velocity moment equations was truncated at second order, and the unclosed collisional source terms were closed using an isotropic Gaussian (Maxwellian) distribution (Levermore & Morokoff 1996;Vié, Doisneau & Massot 2015). Then, by employing Sturm's theorem (Sturm 1829), it was demonstrated that the resulting two-fluid model is hyperbolic for physically realistic values of the model parameters. In comparison to other two-fluid models, novel contributions to the pressure tensor and energy flux (which appear in closed form) arise and play a key role in ensuring hyperbolicity when fluid and particle material densities satisfy ρ f ρ p . Here, we employ the same model formulation, extended to account for the added mass from particle wakes and pseudo-turbulence, to compressible fluid-particle flows with a slip velocity due, e.g. to buoyancy.
Our treatment of added mass is similar to Cook & Harlow (1984) (see appendix A for more details), but generalized to a compressible fluid and a non-constant added-mass function. The latter is required to handle flows wherein the particle-phase volume fraction varies significantly. In our model and in the model of Cook & Harlow (1984), mathematical objectivity is ensured, unlike in other formulations (e.g. Drew, Cheng & Lahey 1979;Massoudi 2002). In the context of kinetic theory, the approach of Cook & Harlow (1984) where the added mass moves with the particle velocity allows us to simply redefine the particle properties without changing the basic form of the kinetic equation governing the velocity distribution function (Fox 2019). Nonetheless, because the fluid in the particle wake is not fixed, but exchanges with the bulk fluid, mass transfer must be included in the kinetic equation to model the convective mass-transfer process. Here, a simple model is employed that depends on a mass-exchange function S a . (See figure 1 for details.) Because the mass-transfer model involves neither spatial nor temporal derivatives, its form does not affect the hyperbolicity of the two-fluid model.
From a kinetic-theory perspective, the added mass of fluid on a particle can be accounted for by defining a particle's volume and mass to include the fluid moving with the particle u f S a u p S a V p V p V a FIGURE 1. Schematic of a particle with its added volume of fluid (i.e. the wake of the particle). The fluid in the wake exchanges mass with the external fluid at a net rate determined by S a . The total particle volume, moving with velocity u p , is V p with sub-volume V p having material density ρ p and added volume V a having material density ρ f . The external fluid with material density ρ f and moving with velocity u f , has volume V f = V − V p . The mass of the particle is m p = ρ p V p + ρ f V a . In terms of the volume fractions, m p = (ρ p α p + ρ f α a )V = ρ e α p V where ρ e is the effective density of the particle with its added mass and α p = α p + α a . Thus, the added volume of fluid moving with the particle velocity is α a V, and the added mass is ρ f α a V. The added-volume fraction must satisfy 0 ≤ α a ≤ α f so it is convenient to define an added-mass function c m by α a = c m α p α f . As the added volume is usually associated with particle wakes, c m can depend on the particle Reynolds number Re p , the particle-phase volume fraction, and other dimensionless parameters needed to describe the flow. In the limit α p → 1, all of the fluid can be assumed to move with the particle so that c m → 1; however, this is not required for hyperbolicity. (Marchisio & Fox 2013), i.e. the fluid in the particle wake (Moore & Balachandar 2019). The total particle mass m p is then employed in the kinetic-theory expressions for the velocity moments. This procedure introduces two volume fractions, namely α p and α p = α p + α a . The former is the usual volume fraction of the particle phase, while the latter includes the volume of the fluid moving with the particles. Naturally, α p ≥ α p and the corresponding fluid-phase volume fractions are α f = 1 − α p and α f = 1 − α p , respectively. A similar decomposition of the fluid-phase variables is used by Osnes et al. (2019) to define a modified slip velocity (u free = α f u fp /α f ) in supersonic gas-particle flows. Using arguments similar to those of Risso (2018) for bubbly flows, these authors also reported that the pseudo-turbulence in the streamwise direction is well approximated by α a u 2 fp /α f , and for fixed α p show that α a decreases with increasing Re p (Osnes et al. 2020).
For the analysis of hyperbolicity, it is convenient to introduce an added-mass function c m defined such that α a = c m α p α f . In principle, c m can be a function of the slip velocity between the two phases (i.e. of the particle Reynolds number Re p = ρ f d p u fp /μ f where d p is the particle diameter and μ f is fluid viscosity), the density ratio Z = ρ p /ρ f , and the volume fraction α p (Zuber 1964;Sangani, Zhang & Prosperetti 1991;Zhang & Prosperetti 1994). However, in the dilute limit, Odar & Hamilton (1964) found experimentally that c m depends only on the acceleration number Ac = u 2 fp /(ad p ) where a is the magnitude of the particle acceleration, which we approximation below using the drag force. Unless c m = 0, the phase velocities found from the kinetic-theory derivation will not be equal to those found from volume averaging unless added mass is accounted for in the latter. Nevertheless, the kinetic-theory derivation leads to conservation laws in the form of hyperbolic equations. This has advantages over a formulation where the virtual mass is treated as an interphase force when solving the two-fluid model numerically. Fox, F. Laurent and A. Vié Finally, because the added mass can vary from location to location in the flow, mass transfer between the bulk fluid and the added-mass fluid (i.e. the particle wake) must be accounted for in the model. This is done by introducing an added-mass exchange rate S a . The exchange of mass between the bulk fluid and added mass also induces an exchange of momentum and kinetic energy, which depends on the direction of the mass exchange. The bulk-fluid momentum is ρ f α f u f , while that of the added-mass fluid is ρ f α a u p . Concerning the total energy, for the particle phase it is defined by where γ p = 5/3 for hard spheres, Θ p is the granular temperature and ρ e α p = ρ p α p + ρ f α a defines ρ e . For simplicity, in (1.1) the internal energy associated with the solid phase and the added mass is neglected. Otherwise, an additional scalar transport equation would be required, which does not change the hyperbolicity of the system (Houim & Oran 2016). For the bulk fluid, the total energy is defined by where γ f is the fluid specific heat ratio, Θ f is the fluid temperature and k f is the pseudo-turbulent kinetic energy (PTKE). In the two-fluid model, the momentum-exchange contribution is equal to S a u f or S a u p , and the energy-exchange contribution to S a (u 2 f /2 + k f ) or S a E p , depending on the sign of S a . The asymmetry in the energy exchange from the bulk fluid to the particle wake results from neglecting the internal energy in (1.1). In the compressible two-fluid model, (1.1) and (1.2) are used to find the temperatures Θ p and Θ f given the total energies E p and E f , respectively. In the stiffened-gas model used for the fluid, Θ f must be initialized such that the fluid pressure p f is positive.
2. Two-fluid model for compressible flows 2.1. Governing equations The governing equations for mono-disperse particles in a compressible fluid with added mass, but neglecting PTKE, are given in table 1. If PTKE is taken into account (Shallcross, Fox & Capecelatro 2020), the model has the form given in table 2. We should point out that in the balance equation for k f the part of the source term D PT due to drag is Ku 2 fp , which is the same as the correlated part of the source term for total energy D E . Physically, this implies that viscous losses are ignored during the exchange process such that drag transfers energy to k f from the particle phase, which is subsequently dissipated to uncorrelated energy by ε PT . The accuracy of this assumption is likely to depend on the particle Reynolds number, i.e. it will be more accurate for high Re p where the particle wakes are turbulent. In practice, this difference can be accounted for by multiplying Ku 2 fp in D PT (but not in D E ) by a damping factor dependent on Re p . Doing so, it may be possible to reduce the Mach number dependence of C f observed in Shallcross et al. (2020).
In prior work (Fox 2019), it has been demonstrated that for an ideal gas (γ f = 5/3) with material densities such that ρ p ρ f and α a = 0 the two-fluid model in table 1 is hyperbolic for physically relevant values of the parameters. In this work, we mainly consider the opposite case with ρ p ρ f where the fluid phase is described by the stiffened-gas model (Harlow & Amsden 1971;Saurel & Abgrall 1999). For a pure fluid, the latter gives an equation of state of the form and the added-mass source terms are The other variables are defined as follows: Compressible two-fluid model for particles in a fluid modelled as a stiffened gas. Typical values of the specific heat ratios are γ f = 29/4 and γ p = 5/3, and for the stiffened-gas constant p o f = 10 8 kg m −1 s 2 : C D is the drag coefficient that depends on the particle Reynolds number Re p , fluid Mach number and volume fraction; and g is gravity. The default added-mass function is c m = min(1 + 2α p , 2)/2.
In addition to the variables defined in table 1, Compressible two-fluid model for particles in a fluid modelled as a stiffened gas with a transport equation for PTKE k f . The pseudo-turbulence tensor R f arises due to the finite size of the particles and b is the PTKE anisotropy tensor (Tenneti, Garg & Subramaniam 2011). The model for a is based on the asymptotic behaviours for ρ f = 0 and ρ p = 0. The parameter b fixes the ratio 3Θ p /2k f when ρ p = 0, and direct numerical simulation data (Tavanashad et al. 2019) suggest that b = 0.365. The constant C f is order one and fixes the magnitude of k f in spatially homogeneous flow (Shallcross et al. 2020). An alternative is to use a transport equation for ε PT to account for the integral length scale of PTKE in lieu of d p .
used to set the speed of sound in the fluid phase. For example, water can be simulated with p o f ≈ 2225 MPa. The fluid temperature Θ f (m 2 s −2 ) is found from the fluid energy E f as shown in table 1, and must be large enough that p f > 0. In this work, we will use a stiffened-gas model of the form The actual value of p o f is not important as long as the speed of sound is much larger than the other characteristic velocities (or eigenvalues) of the system. The factor α f /α f has been added to handle the limiting case where α f → 0 (i.e. densely packed particles), for which this ratio diverges. Other forms of the stiffened-gas model are possible, and the factor is not needed for more dilute systems where the disperse-phase eigenvalues remain well separated from those of the fluid phase. For the disperse (i.e. particle) phase, the radial distribution function g 0 controls the speed of sound as α f approaches zero. For example, if g 0 is replaced with unity, the particle-phase speed of sound is weakly dependent on α p . Here, to analysis the hyperbolicity of the two-fluid model, we use a form for g 0 applicable to non-deformable spheres, but other forms can be used as long as 1 ≤ g 0 . Furthermore, replacing α f /α f with g 0 in (2.1) will not change the conclusions drawn in § 3 concerning the hyperbolicity of the two-fluid models.
The kinetic-theory model derived in Fox (2019) made specific assumptions concerning the two-particle distribution function that may be inaccurate for non-ideal gases and liquids. Specifically, the terms involving R and r in table 1 are exact for hard-sphere collisions (i.e. γ f = γ p = 5/3), but their definition in a stiffened gas is less obvious (e.g. should they depend on both γ f and γ p ?). Thus, in our analysis of the hyperbolicity of the two-fluid model in § 3 we also consider a simplified version where these terms are neglected in the fluid phase. Nonetheless, because the particle-phase pressure tensor P p includes an added-mass contribution involving R, one can argue that P p has its origins in the kinetic-theory description. In fact, in § 3 we show that the eigenvalues of the one-dimensional (1-D) model are mainly determined by the choice of P p and p f , with R and r in the fluid phase only slightly changing the eigenvalues (while making the hyperbolicity analysis more complicated). Thus, from the standpoint of applications to real systems, the simplified model may offer a good compromise between computation cost and model accuracy. However, one would also need to account for inelastic collisions, particle-phase viscosity, as well as other effects (see, e.g. Abbas et al. 2010) in most applications, none of which affect the hyperbolicity.

Added-mass model
In addition to the fluid drag with coefficient K, the models in tables 1 and 2 include the buoyancy force, compressibility, lift and added mass. Compressibility and lift are contained in the exchange force F pf (Fox 2019). The added-mass contribution is treated differently than in most other two-fluid models where balance equations are written for each phase with a virtual-mass force. Instead, here the phases are defined by their velocities u p and u f , and the added mass moves with the particle velocity u p (see discussion in Cook & Harlow 1984). For example, the mass per unit volume of the fluid phase moving with velocity so that the mixture density is independent of α a . The various volume fractions appearing in the model are related by Given the conserved variables (X 1 , X 2 , X 3 ) = (ρ p α p , ρ e α p , ρ f α f ) and the particle density ρ p , the volume fractions and fluid density are uniquely determined by Hereinafter, we define the variable c m such that α a = c m α f α p , which is a convenient form to enforce the upper limit on α a . Although its definition is not required to analyse the hyperbolicity, the added-mass exchange rate will be approximated by a linear relaxation model Physically, τ a is the time scale characterizing the expansion/contraction/formation of particle wakes. For example, when a particle moves from a region with large α p to one with small α p (i.e. to larger spacing between particles), c m will be smaller than c m . Thus, the wake will grow by entraining fluid with velocity u f and kinetic energy u 2 f /2 + k f . The time scale in (2.6) is meant to estimate this rate of growth and can be further refined using data from particle-resolved direct numerical simulations (PRDNS) (see, e.g. Moore & Balachandar 2019).

Added-mass function
In principle, by formulating a physically accurate function for c m , the two-fluid model will be able to account correctly for unsteady effects. For example, c m might depend on Ac (Odar & Hamilton 1964), making the added mass of the particle larger when the particle acceleration is high. In this work, we are primary interested in the effect of added mass on the hyperbolicity of the two-fluid model. In this context, source terms that do not depend on space or time derivatives (such as S a ) have no influence on the eigenvalues of the flux matrix. Nonetheless, the added-mass function c m (x) must have the properties 0 ≤ α p c m ≤ 1 and c m (0) = C m where C m is the added-mass constant, which is equal to 1/2 for a spherical particle when α p = 0 (Zuber 1964). In addition, one might require c m (1) = 1 to force all of the fluid phase to be treated as added mass when its volume fraction approaches zero, but this is not required for hyperbolicity.
Theoretical expressions for the dependence of added mass on the particle volume fraction can be found in Sangani et al. (1991); however, these expressions are valid for α p < 0.5. From the hyperbolicity analysis in § 3, we find that 0.085 < c m < 1/α p , which corresponds physically to 0 < α f < α f . These observations suggest the use of the expression proposed by Zuber (1964) (written to account for the difference between u f and v f ) of the form c m = 1 2 min 1 + 2α p , 2 .
(2.7) Sangani et al. (1991) showed that this form is suitable for most applications (e.g., bubble and spherical particles with no-slip and free-slip boundaries); hence, it will be The diagonal line corresponds to c m = 0. For the function in (2.7), the dependence will be the same as C m = 1 for 1/2 < α p . Note that the curve for C m = 1 is the same for both choices of x.
the default expression in the proposed two-fluid model. Nonetheless, as done in Moore & Balachandar (2019) for the velocity wakes around Lagrangian particles, PRDNS could be used to improve this model to account for the particle Reynolds number and volume fraction dependencies. Another possible expression (which allows for direct computation of α p versus α p ) is the linear form Based on their experimental results, Odar & Hamilton (1964) found that C m depends on the acceleration number as where β ≈ 3 and the acceleration number is defined by Ac = 4/(3C D ). Thus, for very slow acceleration, C m = 1/2, whereas for rapid acceleration C m = 1. However, more recent theoretical works (e.g. Auton, Hunt & Prud'homme 1988;Sangani et al. 1991;Mei & Adrian 1992) suggest that the decomposition of the virtual-mass and history forces used by Odar & Hamilton (1964) is not reliable, and that C m = 1/2 independent of Ac. In any case, using (2.8) and given that (2.10) For C m = 1, the desired root is α 2 f = α f (which also holds for (2.7) when α f < 1/2). For other values of C m , the root can be found numerically as illustrated in figure 2.
As previously noted, the choice of c m has no effect on the hyperbolicity of the two-fluid model. Notwithstanding this fact, for actual applications, it will be important to choose a functional form that accurately matches the dependence of the added mass on α p , etc., derived from PRDNS, experiments and theory.

Particle-phase pressure tensor
In fluid-particle flows, the particles have uncorrelated motion due to fluid-mediated interactions and direct collisions (Lhuillier et al. 2013). In recent PRDNS studies of bubbly flow (du Cluzeau, Bois & Toutant 2019; du Cluzeau et al. 2020), these fluctuations are referred to as the dispersed-phase Reynolds stresses, but it is important to note that they are present in purely laminar flows (Biesheuvel & van Wijngaarden 1984). Indeed, in kinetic theory the magnitude of the dispersed-phase Reynolds stresses is proportional to the granular temperature Θ p . In du Cluzeau et al. (2020), it is found that these terms (referred to as M extra and M LD ) make a significant contribution to the dispersed-phase momentum balance. As described by these authors, in two-fluid models the corresponding flux terms in the dispersed-phase momentum balance are usually separated into 'dispersion' forces (proportional to ∂ x α p ) and a 'drag' force contributions. However, from the standpoint of examining the hyperbolicity of the two-fluid model, it is simplest to treat them as part of the momentum flux as done in this work.
Considering the effective repulsive force exerted between particles in random motion, Batchelor (1988) proposed a (1-D) particle-phase stress model written in terms of the hydrodynamic diffusivity D and the bulk mobility B of the form (written in our notation) where m p is the particle mass. He then used physical reasoning to argue that Considering that Batchelor's model was developed for a 1-D flow with an incompressible fluid phase, it is not unreasonable to treat his dispersion term as part of the particle pressure as done in our compressible two-fluid model. From the kinetic-theory derivation (Fox 2019), a dispersion term is also found in F pf , but written in terms of ∂ x ρ f and not ∂ x α f . Thus, this dispersion term would be zero for an incompressible fluid phase. Mathematically, the dispersion term in (2.11) would appear with the opposite sign in the fluid momentum balance (i.e. it would be an interphase force), and the mixture momentum balance would only contain p p . Syamlal (2011) derived a 'buoyant-force' term that extends the fluid pressure in the Archimedes force to include a relative-velocity contribution of the form ρ f α f u fp ⊗ u fp . Neglecting the particle pressure and considering an incompressible fluid, he demonstrates that the two-fluid model for mass and momentum with this additional force is hyperbolic. Comparing his model to the one in table 1 (and ignoring added mass), we observe that the term in the fluid-phase momentum flux involving R (which is exact from kinetic theory (Fox 2019)) is not present, and the particle-phase pressure tensor term ∂ x · (α a P a fp ) is replaced by the buoyant-force term α p ∂ x · (ρ f α f u fp ⊗ u fp ). In § 3.2, we find that the R contribution to the fluid-phase momentum flux is not required for hyperbolicity (and, for constant ρ f , can be combined with the fluid pressure as done in § 5.2). Thus, by rewriting the buoyant-force term in the form it can be interpreted as a combination of a fluid-mediated particle-phase pressure tensor and a dispersion force, albeit with a negative coefficient. As we show in § 3, the fluid-mediated particle-phase pressure is essential for ensuring a hyperbolic system. Zhang, Ma & Rauenzahn (2006) and Zhang (2020) derived a two-fluid formulation from a general kinetic theory accounting for long-range particle-particle interactions. A particle-fluid-particle (PFP) force of the form ∂ x · (α p Σ pfp ), where Σ pfp is the PFP stress, appears in their formulation. For uniform potential flow with constant ρ f , Zhang (2020) finds that where C 1 and C 2 are coefficients that can be determined numerically. Unlike in (2.13), no dispersion term arises in addition to the PFP force, nor is Σ pfp related to the Archimedes force. Nevertheless, the stress tensor in (2.13) is a special case of (2.14) and, hence, it is reasonable to expect that the PFP force would affect favourably the hyperbolicity of the system (which depends on the trace of Σ pfp ∝ ρ f u 2 fp ). In kinetic theory, the dispersed-phase Reynolds stresses and fluid-mediated interactions contribute to the particle-phase pressure tensor. Thus, from a physical-modelling standpoint, an important component in the two-fluid model is the closure for this term which is a particular form of (2.14). This model for P p combines the kinetic-theory dependence on Θ p due to uncorrelated velocity fluctuations and direct collisions when ρ f ρ p (i.e. p p ) with a component to describe the fluid-mediated interactions between particles that are taken to be proportional to the added mass. In other words, even when the granular temperature is null, in order to have a globally hyperbolic system we assume that a particle pressure exists due to interactions between the particles via the fluid phase (see van Wijngaarden 1976;Batchelor 1988;Zhang et al. 2006;Zhang 2020, for detailed discussions).
In (2.16), the contribution 1 + 4α p g 0 with 1 ≤ g 0 accounts for the excluded volume occupied by the particles. Other formulations of (2.16) are possible and will perhaps be required to capture the correct physics (e.g. when ρ f ≈ ρ p or for deformable particles such as bubbles). For example, one might consider replacing α p g 0 with α p g 0 , or changing altogether the definition of g 0 . However, as shown in § 3, such changes will not affect the hyperbolicity of the compressible two-fluid model as long as γ p is not so large as to make P a fp negligible. Furthermore, in § 3 we find that with α p = 0, the system is hyperbolic even when c m = 0 (i.e. α a = 0 in (2.15)). Thus, it is possible for P a fp in (2.16) to depend linearly on α p (so that the fluid-mediated pressure depends on α 2 p ) without changing the hyperbolicity of the system. An example of this behaviour is presented in appendix B where it is shown that for an incompressible fluid the two-fluid model with trace(α a P a fp ) ∝ (α p ) 2 u 2 fp is globally hyperbolic. The tensorial form of the fluid-mediated particle pressure in (2.16), and its appearance with the opposite sign in the fluid-phase momentum balance, is motivated as follows.
In the kinetic-theory derivation of Fox (2019), it was shown that the mixture pressure tensor has the form P mix = P 1 + P 2 + c 12 P BE regardless of the size ratio between the hard spheres. The Boltzmann-Enskog contribution P BE leads to the term involving R in the fluid-phase momentum balance when added mass is neglected (α a = 0). Biesheuvel & van Wijngaarden (1984) derive a contribution to the mixture stress and liquid-phase Reynolds stresses with the same tensorial form as R based on potential flow around spherical bubbles, but neglecting particle-particle interactions. Thus, with added mass, we assume that the mixture pressure tensor remains unchanged, and share the contribution α a P a fp between the two phases. This is consistent with the kinetic-theory derivation where the particle pressures in each phase depend on the particle size ratio, while the mixture pressure does not (Fox 2019). Finally, note that the contribution ∂ x · (α a P a fp ) arises from the kinetic-theory derivation as a modification to the pressure tensors, while in the compressible two-fluid model it can be interpreted as a fluid-mediated exchange force. Finally, the parameter γ p in (2.16) is equal to 5/3 for a hard sphere in an ideal gas, but in general it can be used as a parameter to set the magnitude of the fluid-mediated particle pressure (i.e. tr(P a fp ) ∝ 5/(3γ p )). On the other hand, the tensorial form of P a fp must be kept consistent with that of R as both arise due from the same term in the kinetic-theory derivation (Fox 2019). As seen from (2.14), up to the scalar coefficients that can depend on α p , this tensorial form is the only one that can be formed from u fp (Zhang 2020).
In summary, the particle-pressure tensor in (2.15) combines two limiting behaviours for the material-density ratio and it is a key modelling component for ensuring hyperbolicity when ρ p ρ f . This is consistent with Batchelor (1988) where it is also shown to have a strong effect on the linear stability of a uniform suspension. It is also consistent with the kinetic-theory derivation of Zhang et al. (2006); Zhang (2020) who demonstrate that an inter-species stress arises due to particle-fluid-particle interactions, which do not depend on Θ p . Nonetheless, future research should be devoted to refining the model for α a P a fp in (2.16) to account for the α p and Re p dependencies observed in PRDNS.

Limiting cases
Having previously investigated the case with ρ f ρ p (Fox 2019), in the remainder of this work we are particularly interested in the limiting cases with ρ p = 0 (i.e. the particles have zero mass) so that ρ e α p = ρ f α a . However, we will also analyse the hyperbolicity of the complete model for selected values of the material density ratio Z. As is well known, the drag and body forces appearing on the right-hand sides of the model equations do not affect the eigenvalues of the two-fluid model. Hence, in § 3 we will ignore them and consider only the terms involving temporal and spatial gradients.

Hyperbolicity of 1-D two-fluid model
In order to determine whether the full three-dimensional (3-D) model in table 2 is hyperbolic, it suffices to consider a system with one spatial direction (see, e.g. Ndjinga 2007; Kumbaro & Ndjinga 2011;Lhuillier et al. 2013). This approach will be followed here, starting with the 1-D model without source terms that do not involve temporal or spatial derivatives.

One-dimensional compressible two-fluid model
The 1-D model without the source terms is given in table 3, written in terms of eight independent variables We define Z such that ρ f = Zρ p . As ρ p is constant, it can be factored out of the model if desired. The conserved variables (X 1 , X 2 , X 3 ) are related to (α p , Z, α f ) by and all other variables appearing in the momentum and energy balances can be found from these variables. In addition to the model in table 3, we will also analyse the simplified model given in table 4, which neglects the Boltzmann-Enskog fluxes (i.e. R and r in the fluid phase) and forces (i.e. F fp and D fp ) that are specific to hard-sphere mixtures (Fox 2019). Formally, the 1-D models can be rewritten as Here, B 0 is the contribution due to the pressure and buoyancy terms, and A and B * are  3. This model is hyperbolic when the fluid-phase eigenvalues are sufficiently separated from those for the particle phase. When this is not the case, the kinetic-theory terms in the full model may be needed to achieve hyperbolicity.
from the convection terms. Written in terms of the components of X , the latter are and (3.5) The components of B 0 are more complex due to the nonlinearities, but can easily be computed using symbolic software, as can the flux matrix and its characteristic polynomial. Due to the nonlinearities of the additional fluxes and forces in the full versus the simplified model, the latter can be analysed analytically in greater detail. Nonetheless, it is always possible to compute the eigenvalues numerically in order to check the predictions of the analysis. The eight eigenvalues of F can be written u f + u 0 λ k with k ∈ {1, . . . , 8}, where for fixed values of p o f = Z 0 p o /u 2 0 = p o f /ρ p u 2 0 , γ f and γ p , each λ k , called here a normalized eigenvalue, depends on five dimensionless parameters where Θ f = u 2 0 + Θ 0 and Θ 0 is defined by p f = 0 from the stiffened-gas model. The parameter Θ 0 depends on Z. The λ k are the roots of P(X) = Q(u 0 X)/u 8 0 , where Q is the characteristic polynomial of F − u f I. In general, in order for the eigenvalues to be real, p f must be positive. The characteristic velocity u 0 should not be confused with the speed of sound in the stiffened-gas model, which scales like c f = (γ f p o f ) 1/2 and is orders of magnitude larger than u 0 . For the model in table 3, there are two normalized eigenvalues that can be computed analytically, namely, 0 and Ma s . For the model in table 4, there is an additional normalized eigenvalue at Ma s . In general, the remaining normalized eigenvalues in both models depend on the five parameters in (3.6a-e).
For α p = 0, the normalized eigenvalues (which are the same for both models) can be computed analytically and are always real valued, including when c m = 0. Here, the two 'fluid-phase' eigenvalues that depend on p o f are always real and distinct. Note that when Θ r = 0 the 'particle-phase' eigenvalues scale with Ma s , but always remain real-valued. When Ma s = 0, these eigenvalues depend on Θ r like an ideal gas (γ p = 5/3). In the following, we investigate the behaviour of the eigenvalues for a stiffened gas (γ f = 29/4) with fixed values of Z, namely, +∞, 1, and 0; which correspond, respectively, to bubbly, neutrally buoyant and granular flow. The behaviour of the eigenvalues for other values of Z can be inferred from these limiting cases. For the model in table 3, there will be six eigenvalues that vary with α p , as opposed to five for the model in table 4. From the examples in figure 3, it can be observed that the differences between the full and simplified model are small. The magnitudes of the two 'particle-phase' eigenvalues increase with α p mainly due to g 0 , while their values at α p = 0 depend on c m as shown in (3.7a-d).

Hyperbolicity analysis of simplified model
For the simplified model in table 4, a theoretical study of the hyperbolicity can be carried out analytically. As done in Chalons et al. (2017), Sturm's theorem (Sturm 1829) can be used, which determines the number of distinct real roots in a given interval. For that, let us consider the polynomial P 0 = P/(X(X − Ma s ) 2 ), where P is the polynomial defined above. The Sturm sequence of polynomials is defined by P 0 , P 1 = P 0 and, for any n ∈ {0, 1, 2, 3}, −P n+2 is the remainder of the Euclidean division of P n+1 by P n . With the use of symbolic software, one can compute this sequence. If the coefficient S n of the highest-order term of each P n , called hereinafter a Sturm's coefficient, is positive for n ∈ {0, 1, . . . , 5}, then Q has five real roots, meaning that all eigenvalues of the system are real. In the general case, it is hard to prove that all the Sturm's coefficients S n are positive. However, since c f is large compared to u 0 , the limit when c f tends to infinity can be studied, i.e. for a very large value of p o f . Thus, a Taylor expansion of the S n can be done when = 1/ p o f tends to zero, for γ f = 29/4 and γ p = 5/3. Then, S 0 = 1, S 1 = 6 and the limit of S n when tends to zero is studied 10)

Eigenvalues for specific cases
We are specifically interested to know whether the 1-D models are hyperbolic for all physically relevant values of 0 ≤ α p ≤ 1 and 0 ≤ c m ≤ 1/α p . As can be seen in (3.7a-d), K r mainly affects the 'fluid-phase' eigenvalues, so it suffices to show hyperbolicity for K r = 0. Similarly, Θ r mainly affects the 'particle-phase' eigenvalues and Θ r = 0 is known to yield complex eigenvalues when added mass is neglected (Fox 2019); hence, the analysis of this limiting case is of particular interest. As done in Fox (2019), we will make use of stability plots found from the Sturm coefficients to check for complex eigenvalues in α p -c m parameter space.
3.3.1. Granular flow with ρ f = 0 The limit Z = 0 corresponds to a granular flow with ρ f = 0. For this case, the 1-D model is globally hyperbolic. Nonetheless, the hyperbolicity plot in figure 4 requires a Ma s -dependent minimum value for Θ r due to round-off errors in evaluating the Sturm coefficients. As can be seen in figure 3 and from (3.7a-d), the particle phase has multiple eigenvalues at Ma s when Z = 0 and Θ r = 0.

Neutrally buoyant flow with ρ p = ρ f
The limit Z = 1 corresponds to a neutrally buoyant flow with ρ p = ρ f . From the hyperbolicity plot in figure 5, we can observe that the models are hyperbolic except for a small region near c m = 0. In other words, there is a minimum value of c m above which the models are globally hyperbolic. Note that this value is significantly smaller than the standard added-mass constant c m = 1/2. In figure 6, examples with complex eigenvalues corresponding to small c m are shown. It can be observed that with Z = 1 the eigenvalues for the two models are quite similar unless c m is very small. It is noteworthy that when c m is small enough, the 'disturbance' eigenvalue that begins above unity at α p = 0 can be less that unity for values of α p near 0.15. In other words, particle-phase disturbances propagate more slowly than the mean slip velocity if the added mass is relatively small.

3.3.3.
Bubbly flow with ρ p = 0 The limit Z → ∞ corresponds to a bubbly flow with ρ p = 0. From the hyperbolicity plot in figure 7, we can again observe that the models are hyperbolic except for a small region near c m = 0. Furthermore, for the simplified model, the non-hyperbolic region is very small and can be easily avoided by proper choices for c m and τ a . In figure 6, examples with   complex eigenvalues corresponding to small c m are shown. It can be observed that with Z = 10 000 the eigenvalues for the two models are nearly identical. In general, as predicted from the hyperbolicity plot in figure 7, the full model has the largest region of parameter space in which it is hyperbolic. As mentioned earlier, the Boltzmann-Enskog fluxes are valid for a hard-sphere mixture, so neglecting them as done in the simplified model may be  allowable without significantly changing the hyperbolicity. The added-mass contribution to the particle-phase pressure tensor then determines the domain of hyperbolicity of the two-fluid model. In conclusion, it is worth noting that in the full model the region of hyperbolicity for 0.7 < α p includes c m → 0. In other words, when the particle-phase volume fraction is near unity, the added mass has no effect on the well posedness of the full model. Thus, c m can take any value in the interval [0, 1] for densely packed particles. In general, the effect of added mass on hyperbolicity is most important for α p ≈ 0.1, regardless of the material-density ratio.

Numerical examples of 1-D model
To illustrate the behaviour of the proposed two-fluid model, in this section we develop a 1-D numerical solver for the full model written in conservative form. In table 5, the 1-D

Conservative form of 1-D model
The added-mass contribution to the particle-phase pressure tensor P a fp is purely mechanical. This implies that the granular energy balance for Θ p has a compression source term that depends only on the 'thermodynamic' pressure p p (see appendix A in Houim & Oran 2016, for details). Without this property, the source term can generate a negative granular temperature when the thermodynamic pressure is null. Using the 1-D model in table 5, the granular energy balance becomes (with γ p = 5/3) which has the necessary property that the right-hand side is non-negative when Θ p is null. In contrast, if p p were replaced by ( p p + P a fp ) in (4.1), then the first term on the right-hand side would be negative when Θ p = 0 and ∂ x X 4 is positive (i.e. during expansion of the particle phase), leading to a non-physical negative granular temperature. Finally, note that body forces (i.e. gravity) do not appear in (4.1) and, therefore, it can be solved in place of the total energy balance to avoid the associated numerical errors observed in § 4.3.
The other model parameters are as follows: 3/2 8 c m = 1 2 min 1 + 2α p , 2 TABLE 5. One-dimensional compressible two-fluid model in conservative form with densities and pressures normalized by the particle material density ρ p . The terms on the left-hand side are the conservative fluxes, while those on the right are interphase exchange terms and g x is the component of gravity in the x direction. The fluid kinematic viscosity is ν f and the particle diameter is d p . For water, ν f = 10 −6 m 2 s −1 and the stiffened-gas model parameters are γ f = 29/4 and p o = 10 8 m 2 s −2 . Z 0 is the reference density ratio and C f = 1. For the particle phase, γ p = 5/3 and C D is the Re p -dependent drag coefficient where, for Stokes drag, C D Re p = 24.
Hyperbolic compressible two-fluid model 903 A5-23 The mixture mass ( = Y 2 + Y 3 ), momentum M = Y 4 + Y 5 and energy E = Y 6 + Y 7 balances can be written as which have the form of hyperbolic conservation laws, albeit with rather complex momentum and energy fluxes. From a computational standpoint, (4.2) can be solved with the PTKE and particle-phase balances using operator splitting for the left-/right-hand sides, respectively. Here, the left-hand sides are fluxes in the finite-volume sense, whereas the right-hand sides are source terms. Alternatively, the balance for Y 6 can be replaced with the granular energy balance for Θ p shown in (4.1), which is preferable for dissipative systems to avoid round-off errors due to the very small value of Θ p (see Houim & Oran (2016), for a discussion of this point), and for the numerical treatment of body forces.

Numerical solver
The 1-D model in table 5 has the form where f are the spatial fluxes, and h are the source terms. In the numerical solver for (4.4), the conservative variables in each grid cell are updated using an explicit algorithm The source term in (4.5) is evaluated using a central-difference formula for the spatial gradient. In principle, the source term could be stiff and one might want to use an implicit solver. However, we found that the time step required for the spatial fluxes is small enough that this is unnecessary. We should note that for simulations starting with zero particle-phase energy, the explicit Euler scheme for the gravity term generates a negative granular temperature (i.e. Y 0 4 = Y 0 6 = 0 yields Y 1 4 = tY 0 2 g x , Y 1 6 = 0). In general, if there is mean slip between the phases (i.e. Z 0 / = 1), the granular temperature becomes positive everywhere in the domain after a few time steps. To avoid such numerical issues, when evaluating the particle pressure and spatial fluxes, the granular temperature is set to zero whenever it is negative.

A5-24
R. O. Fox, F. Laurent and A. Vié The numerical fluxes in (4.5) are defined using a classical HLL approach (Harten, Lax & van Leer 1983;Toro 1997) where a − < 0 < a + correspond to the minimum and maximum eigenvalues of the system. In all cases considered below, these eigenvalues come from the stiffened-gas model and a − ≈ −a + . In our simulations, the eigenvalues are computed at each time step from the characteristic polynomial. On a uniform grid with spacing x, the time step is set using x/ t = 2 max(a + , −a − ). As described in § 3, the six other eigenvalues of the 1-D system are order one in magnitude, and thus are at least two orders of magnitude smaller than a + . As a consequence, the HLL fluxes in (4.6), designed for two-wave systems (Toro 1997), will generate significant numerical diffusion for the system in table 5. For example, the effective diffusivity of the volume fraction is D ∝ x a + due to the final difference term on the right-hand side of (4.6). For this reason, unless x is very small, the material interface for the density-matched case (Z 0 = 1), for which the mean velocity is very small, will be smeared out over time. Future work should therefore focus on developing hyperbolic solvers specifically for multiphase systems to reduce the numerical diffusion using higher-order spatial reconstruction schemes (Toro 1997).

Numerical examples
The numerical examples provided in this section illustrate the behaviour of the model for Riemann problems with different initial conditions on the right/left sides of the domain. In all cases, the initial value Z = Z 0 is used (i.e. fixed material-density ratio) and corresponds to monodisperse particles in a given fluid with kinematic viscosity ν f = 10 −6 m 2 s −1 . The fluid temperature Θ f is set to be 5000 m 2 s −2 larger than Θ 0 in the stiffened-gas model. For simplicity, we consider Stokes drag (C D Re p = 24), a particle diameter of d p = 10 −3 m and set c m = 1/2. Finally, in order to keep the time step reasonable, we set p o = 10 4 m 2 s −2 , which yields a + ≈ 775 m s −1 . Increasing p o will result in smaller variations in Z, but requires a correspondingly smaller time step. The computational domain is taken as x ∈ [−1/2, 1/2] m and zero-flux boundary conditions are employed on each end. To illustrate the effect of numerical diffusion in the HLL scheme, the grid spacing is set to x = 1/N m with N = (1000, 2000, 4000).
In the first example, we consider a case with Z 0 = 1. The initial conditions are α p = 0 on the left half and α p = 0.1 on the right half of the domain. The fluid and particle velocities and the granular temperature are null. With gravity, a pressure gradient develops in the fluid phase and Z becomes very weakly dependent on x. Because the particles have the same density as the fluid, the exact solution for α p does not change with time. However, as seen in figure 8, the HLL fluxes are diffusive, leading to a smearing out of the material interface. Without gravity, numerical diffusion is also observed for α p even though the velocities are null. At shorter times (but still long compared to the fluid speed of sound), the numerical diffusion is less obvious. As expected for Z 0 = 1, aside from the volume fractions and P f , the primitive variables are nearly uniform. Note that the particle pressure P p is very small due to the negligible slip velocity and small Θ p . For the exact solution, Θ p is null and, as discussed earlier, negative values are computed due to the treatment of body forces with the Euler time step. Finally, the added-mass c m remains close to the equilibrium value of 0.5, and thus well above the minimum value required for hyperbolicity.
In the second example, we consider a case with Z 0 = 10 4 corresponding to buoyant particles. The initial conditions are again α p = 0 on the left half and α p = 0.1 on the right  half of the domain. The fluid and particle velocities and the granular temperature are null. With gravity, a pressure gradient develops in the fluid phase and Z becomes very weakly dependent on x. The buoyancy force makes u p positive, moving particles towards the top of the domain. However, as seen in figure 9, the HLL fluxes lead to a smearing out of the material interface, with the numerical diffusion resisting the rise velocity. A finer grid exhibits less numerical diffusion, but the smearing of the volume fraction is still obvious. Presumably, if the grid were made fine enough, α p would remain near zero for x < 0, and be larger at x = 0.5 due to buoyancy. In general, the primitive variables are non-uniform due to the buoyancy force. The particle pressure P p is significant due to the slip velocity. As in the first example, the added-mass c m remains close to the equilibrium value of 0.5, and thus well above the minimum value required for hyperbolicity. A case with Z 0 = 10 −4 (not shown) exhibits similar behaviour, but with the particle volume fraction larger at the bottom of the domain. In summary, aside from excessive numerical diffusion due to the HLL fluxes, solutions to the model in table 5 exhibit the expected physical behaviour. Depending on the value of Z 0 , the particle pressure can play a significant role in the momentum balances. In real applications, Θ p (and, hence, p p ) will be much smaller due to inelastic collisions and lubrication effects (Abbas et al. 2010). However, this will not affect the hyperbolicity. Furthermore, the added-mass c m can vary spatially due to transport between regions with different volume fractions, but always remains well above the minimum value of 0.085 needed to ensure hyperbolicity.
. Primitive variables at t = 0.1 s with Z 0 = 10 4 and x = 1/N m (blue, N = 1000; red, N = 2000; gold, N = 4000). Due to the finer mesh, numerical diffusion is less important for larger N as can be seen from the spatial distribution of α p .

Discussion and conclusions
The compressible two-fluid model in table 2 describes inviscid fluids with arbitrary material-density ratios. As seen from the examples in the previous sections, the model for P p plays a key role in the hyperbolicity of the two-fluid model. In particular, when the material-density ratio Z is non-zero, P p must be non-zero when Θ p = 0 in order to eliminate complex eigenvalues. Following Batchelor (1988) and Zhang et al. (2006), we have thus included an added-mass contribution to the particle-phase pressure tensor that depends on the slip velocity between the phases. For bubbly flow where the particle shape is flexible, the g 0 expression for solid particles is likely to diverge too quickly with increasing α p . On the other hand, with rigid particles a frictional component must be added to P p , which is independent of Θ p . Houim & Oran (2016) have analysed such a case and showed that the eigenvalues remain real valued. Therefore, the overall conclusion is that the two-fluid model in table 2 provides a hyperbolic inviscid model for describing compressible disperse-phase flows for all material-density ratios.

Extension to viscous flows and other interphase forces
The model equations in table 2 can be augmented in different directions. First, to model viscous flows (Abbas et al. 2010;Guazzelli & Pouliquen 2018), (traceless) viscous-stress tensors for the fluid and particle phases can be added to the momentum and energy balances. Note that for the fluid phase, b acts like a pseudo-turbulent viscous stress and can be modelled as a Newtonian fluid with an effective viscosity depending on k f and ε PT . The parameter a appearing in D E and D PT determines the distribution of pseudo-turbulent kinetic equation between k f and Θ p . The latter has been investigated for spatially homogeneous, incompressible flow by Tavanashad et al. (2019) over a wide range of material-density ratios, and these results could be use to develop a correlation for a. Likewise, C f fixes the value of k f /u 2 fp and can be fit to the data of Tavanashad et al. (2019). As with all two-fluid models, a closure for the drag coefficient K must be provided, which will depend, as usual, on the particle Reynolds number and volume fraction in addition to the fluid Mach number. Finally, additional interphase forces can be added to the momentum and energy balances to describe the effects of mean shear and vorticity on the disperse phase. As mentioned earlier, although these forces contain spatial gradients of the phase velocities, they act normal to the flow direction and, therefore, do not change the hyperbolicity of the system. Note that the effect of 'turbulent dispersion' of the disperse phase is already included using the tensor R f , and thus no additional terms are needed in the balances in table 2. The same is true for the virtual-mass force, which is accounted for as added mass. The coefficient of the lift force in F fp should be revisited to account for surface-tension effects with deformable particles (du Cluzeau et al. 2020).

Two-fluid model for bubbly flow with constant ρ f
For bubbly flow where ρ p ρ f , the two-fluid model in table 2 can be simplified by removing the transport equation for E f and treating ρ f as constant. The fluid pressure is then found using the condition that α f + α p = 1. This seven-equation model is shown in table 6. It is important to note that while the mean velocity fields are weakly coupled with the pseudo-turbulence kinetic energy, it is often necessary to solve for k f and Θ p for other purposes. For example, k f will be needed to model the effective viscosity or the effective diffusivity of a passive scalar transported by the fluid (Peng et al. 2019). As mentioned above, the bubbly flow model in table 6 can be augmented with a viscous-stress tensor (including pseudo-turbulence) for the fluid phase. Unlike in most other hyperbolic formulations for bubbly flows (see, for example, Panicker et al. 2018), it is not necessary to add a turbulent-dispersion term to enforce hyperbolicity. As shown in appendix B, the model for P p ensures global hyperbolicity. In the terminology of two-fluid models (see Lhuillier et al. 2013), the seven-equation two-fluid model in table 6 is a two-pressure model with mixture pressure tensor P = ( p p + p f )I. As we show in appendix B, the shared-pressure model with P p = 0 is not hyperbolic (Drew & Passman 1998) and, thus, cannot be used for bubbly flow simulations because it produces non-physical solutions (see examples in Panicker et al. 2018). Lhuillier et al. (2013) discuss the history of effective-field models for disperse flows, providing insights into why past formulations are mathematically ill posed. It is therefore of interest to compare the two-fluid model in table 2 to their formulation in order to understand why it is hyperbolic. However, even before performing this exercise, it is noteworthy that these authors suggest that 'a promising direction is to associate added-mass and the pseudo-turbulence of the particles'. For clarity, in this section we will use the notation developed in this work. However, we should emphasize that as discussed in appendix A the effective-field model is written in terms of the velocity v k , while here we use u k to account for the added mass.

Relation to effective-field models
The pseudo-turbulent kinetic energy in Lhuillier et al. (2013) is denoted by K k , so that K f = k f and K p = (1/(γ p − 1))Θ p in our notation. In other words, as the particles have no internal energy, the granular temperature plays the role of the pseudo-turbulent kinetic energy of the particle phase. The momentum balances in Lhuillier et al. (2013) are written TABLE 6. Seven-equation two-fluid model for bubbly flow with constant fluid density and γ p = 5/3. C D is the drag coefficient that depends on the particle Reynolds number and volume fraction, and g is gravity. The energy balance for E p can be rewritten in terms of Θ p . In principle, this model can be applied for any value of Z provided ρ f can be treated as constant (i.e. low Mach number flows). In the fluid momentum balance, a traceless stress tensor due to R and R f can be included without changing the hyperbolicity of the system. in our notation as (5.1a,b) where D k /Dt is the convected derivative with velocity u k , Π k are the phasic stresses and F = Ku pf is the interphase force (i.e. drag). Comparing with table 2, we observe that Π p = P p and Π f = α p ρ f R − α a P a fp , and that the kinetic-theory expression for P f includes the pseudo-turbulent pressure due to R f . However, as shown earlier, the two-fluid model is hyperbolic even with R f = 0 for the fluid phase, so the most important term is the added-mass-dependent contribution to Π p and Π f (which can alternatively be treated as part of F ).
The energy equation for the particle phase found from the balances in table 2 is where the left-hand side is a non-dissipative pseudo-turbulent kinetic energy exchange term. The parameter a determines the amount of mean kinetic energy that is directly dissipated to fluid-phase internal energy, so for a non-dissipative system a = 0. Recalling that b and ε PT arise due to dissipation of fluid-phase pseudo-turbulent kinetic energy to fluid-phase internal energy, the non-dissipative terms in the pseudo-turbulent kinetic energy balance yield Then, as could be anticipated from the fact that α p ρ e E p + α f ρ f E f obeys a conservation equation, the sum of (5.2) and (5.3) satisfies the energy conservation condition (2.5) in Lhuillier et al. (2013). Nonetheless, it is important to point out that the trace of P p is not the inter-facial pressure of the fluid at the particle surface. Indeed, the Θ p -dependent part of P p arises in kinetic theory due to particles having different instantaneous velocities. Thus, at best, only the α a -dependent part of P p might be assigned to the inter-facial pressure (see Batchelor 1988, for a discussion of the physical reasoning on why this is incorrect). Supporting the arguments made by Lhuillier et al. (2013) (and consistent with Batchelor 1988), taken as a whole these observations suggest that the Θ p -independent contribution to P p is a necessary condition for hyperbolicity of two-fluid models.

Closing remarks
Starting from the kinetic-theory-based model derived from first principles in Fox (2019), the definition of the particle mass was extended to include the added mass moving with the velocity of the particle. This resulted in the compressible two-fluid model in table 1. Then, by relaxing the assumption that the pseudo-turbulent kinetic energy in the fluid phase (denoted by R f ) attain instantaneously its steady-state value, a transport equation was introduced to model its trace (2k f ). The resulting compressible two-fluid model, presented in table 2, has governing equations for pseudo-turbulent kinetic energy in both phases, as well as balance equations for the total energies. The fluid phase is treated as compressible with a stiffened-gas equation of state to describe liquids. As written, the compressible two-fluid model is applicable to flows with an arbitrary material-density ratio Z = ρ f /ρ p .
While needed for accurate physical modelling (e.g. in gas-particle flows), from the hyperbolicity analysis of the 1-D model it was found that the pseudo-turbulent kinetic energies play no role in determining whether the two-fluid model is hyperbolic. In contrast, g 0 and the particle-fluid-particle stress contribution (i.e. α a P a fp ) to P p are crucial for obtaining a hyperbolic model for large to intermediate values of Z. Indeed, for ρ p = 0 (mass-less particles), without these contributions the two-fluid model loses hyperbolicity in physically important regions of parameter space (e.g. Θ p near zero). Future work should therefore focus on obtaining a more fundamental understanding of how to model P a fp and g 0 in the particle-phase pressure tensor for real physical systems, especially for Z ≈ 1. To this end, direct numerical simulations of particle suspensions over a wide range of material-density ratios, Reynolds numbers and volume fractions would be useful (such as is done in Moore & Balachandar 2019;Tavanashad et al. 2019;du Cluzeau et al. 2020), especially if one can unequivocally relate model variables such as c m and P a fp to the data from such simulations (Zhang 2020). Finally, work along the lines of Gu et al. (2019) and Abbas et al. (2010) will be required to account for viscous effects in the particle phase. . Eigenvalues of incompressible bubbly flow model versus α p with α a = c m α f α p , Θ r = 2/γ 2 p and γ p = 5/3. For c d = 1 the eigenvalues are found with g 0 , whereas for c d = 0 the eigenvalues are found with g 0 = 0. Unlike in the compressible model, g 0 is not required to keep the system hyperbolic for large α p . This 'equilibrium' value for Θ r results in one eigenvalue at zero when g 0 = 0, which corresponds to the fluid velocity. Conversely, the equilibrium value for Θ r found from direct-numerical simulation could be used to specific γ p for bubbly flow (i.e. γ p = 4.714 using data from Tavanashad  are real valued for all α p with c m = 1/2, including with g 0 = 0. It is noteworthy that the particle-phase eigenvalues from (3.7a-d) in the limit Z → +∞ are equal to (B 7). This would not be the case if the transport equation for E p were replaced by an algebraic expression for Θ p . In any case, it is remarkable that by adding a transport equation for the added mass and a model for the particle-pressure tensor that does not depend on granular temperature, the incompressible two-fluid model becomes globally hyperbolic. Alternative forms for the particle-pressure tensor are also possible. For example, it is not necessary for α a P a fp to depend on g 0 or α a . In figure 11, the eigenvalues for a fluid-mediated particle-pressure model that is quadratic in α p P p = p p + c f γ p (α p ) 2 α f u 2 fp (B 8) are plotted versus α p . For α p = 0, the two non-constant eigenvalues with this particle-pressure model are 1 ± γ p Θ r , (B 9) and they remain real valued for all α p ∈ [0, 1] if 0.1 ≤ c f , independent of c m . Comparing figures 10 and 11, we observe that with 0 < c m the quadratic dependence on α p causes the two eigenvalues to be equal at α p = 0 when Θ r = 0. (Here, c m = 0 is a singular case where the lower eigenvalue jumps to zero as α p increase from zero.) For c f < 0.1, the eigenvalues are complex for small α p , before becoming real valued for larger volume fractions. Replacing α p by α p in (B 8) gives qualitatively equivalent results, but the value of c f must be slightly larger to ensure hyperbolicity. Although it introduces a new parameter c f into the hyperbolicity analysis, the form of (B 8) is physically motivated by the fact that binary interactions between particles mediated by the fluid scale with α 2 p in a dilute system. Also, note that the partial derivative of the second term in (B 8) with respect to α p has the form of the 'turbulent-dispersion' force that is used to make bubbly flow models hyperbolic (see, e.g. Panicker et al. 2018).
Finally, it is noteworthy that the partial derivatives of P p appearing in the fourth row of B in (B 5) could be interpreted as arising due to separate forces. For example, ∂P p /∂X 2 might be attributed to 'turbulent dispersion', while ∂p p /∂X 6 acts like a pseudo-turbulent turbophoresis. Nonetheless, they all have a common origin in the particle-phase pressure tensor.