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 Reference Shallcross, Fox and Capecelatro2020), 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 $K u_{fp}^2$, 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 $\varepsilon _{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 $K u_{fp}^2$ 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. (Reference Shallcross, Fox and Capecelatro2020).

Table 1. Compressible two-fluid model for particles in a fluid modelled as a stiffened gas. Typical values of the specific heat ratios are $\gamma _f = 29/4$ and $\gamma _p = {5}/{3}$, and for the stiffened-gas constant $p^o_f = 10^8\ \textrm {kg}\,\textrm {m}^{-1}\,\textrm {s}^{2}$: $C_D$ is the drag coefficient that depends on the particle Reynolds number $Re_p$, fluid Mach number and volume fraction; and $\boldsymbol {g}$ is gravity. The default added-mass function is $c_m^{\star } = \min (1 + 2 \alpha _p , 2)/2$.

Table 2. 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 $\boldsymbol {R}_f$ arises due to the finite size of the particles and $\boldsymbol {b}$ is the PTKE anisotropy tensor (Tenneti, Garg & Subramaniam Reference Tenneti, Garg and Subramaniam2011). The model for $a$ is based on the asymptotic behaviours for $\rho _f=0$ and $\rho _p=0$. The parameter $b$ fixes the ratio $3 \varTheta _p / 2 k_f$ when $\rho _p=0$, and direct numerical simulation data (Tavanashad et al. Reference Tavanashad, Passalacqua, Fox and Subramaniam2019) 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. Reference Shallcross, Fox and Capecelatro2020). An alternative is to use a transport equation for $\varepsilon _{PT}$ to account for the integral length scale of PTKE in lieu of $d_p$.

In prior work (Fox Reference Fox2019), it has been demonstrated that for an ideal gas ($\gamma _f=5/3$) with material densities such that $\rho _p \gg \rho _f$ and $\alpha _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 $\rho _p \ll \rho _f$ where the fluid phase is described by the stiffened-gas model (Harlow & Amsden Reference Harlow and Amsden1971; Saurel & Abgrall Reference Saurel and Abgrall1999). For a pure fluid, the latter gives an equation of state of the form $p_f = \rho _f \varTheta _f - p_f^o$ where the constant $p_f^o$ is used to set the speed of sound in the fluid phase. For example, water can be simulated with $p_f^o \approx 2225\ \textrm {MPa}$. The fluid temperature $\varTheta _f$$(\textrm {m}^{2}\,\textrm {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

(2.1)\begin{equation} p_f = \rho_f \varTheta_f - \gamma_f (\gamma_f -1) p^o_f \frac{\alpha_f}{\alpha_f^{\star}} . \end{equation}

The actual value of $p_f^o$ 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 $\alpha _f / \alpha _f^{\star }$ has been added to handle the limiting case where $\alpha _f \to 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 $\alpha _f$ approaches zero. For example, if $g_0$ is replaced with unity, the particle-phase speed of sound is weakly dependent on $\alpha _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 \le g_0$. Furthermore, replacing $\alpha _f / \alpha _f^{\star }$ 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 (Reference Fox2019) made specific assumptions concerning the two-particle distribution function that may be inaccurate for non-ideal gases and liquids. Specifically, the terms involving $\boldsymbol {R}$ and $\boldsymbol {r}$ in table 1 are exact for hard-sphere collisions (i.e. $\gamma _f = \gamma _p = 5/3$), but their definition in a stiffened gas is less obvious (e.g. should they depend on both $\gamma _f$ and $\gamma _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 $\boldsymbol {P}_p$ includes an added-mass contribution involving $\boldsymbol {R}$, one can argue that $\boldsymbol {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 $\boldsymbol {P}_p$ and $p_f$, with $\boldsymbol {R}$ and $\boldsymbol {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. Reference Abbas, Climent, Parmentier and Simonin2010) in most applications, none of which affect the hyperbolicity.

2.2. 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 $\boldsymbol {F}_{pf}$ (Fox Reference Fox2019). 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 $\boldsymbol {u}_p$ and $\boldsymbol {u}_f$, and the added mass moves with the particle velocity $\boldsymbol {u}_p$ (see discussion in Cook & Harlow Reference Cook and Harlow1984). For example, the mass per unit volume of the fluid phase moving with velocity $\boldsymbol {u}_f$ is $\rho _f \alpha _f^{\star } = \rho _f ( \alpha _f - \alpha _a )$. Note that

(2.2)\begin{equation} \rho_f \alpha_f^{\star} + \rho_e \alpha_p^{\star} = \rho_f \alpha_f +\rho_p \alpha_p, \end{equation}

so that the mixture density is independent of $\alpha _a$. The various volume fractions appearing in the model are related by

(2.3a–c)\begin{equation} \alpha_f^{\star} = \alpha_f - \alpha_a \quad \alpha_p^{\star} = \alpha_p + \alpha_a \quad \alpha_p + \alpha_f = 1. \end{equation}

Given the conserved variables $(X_1,X_2,X_3)=(\rho _p \alpha _p , \rho _e \alpha _p^{\star }, \rho _f \alpha _f^{\star })$ and the particle density $\rho _p$, the volume fractions and fluid density are uniquely determined by

(2.4)\begin{equation} (\alpha_p, \rho_f, \alpha_a) = \left( \frac{X_1}{\rho_p} , \frac{X_3+X_2-X_1}{\alpha_f}, \frac{X_2 - X_1}{\rho_f} \right) . \end{equation}

Hereinafter, we define the variable $c_m$ such that $\alpha _a = c_m \alpha _f \alpha _p$, which is a convenient form to enforce the upper limit on $\alpha _a$.

Although its definition is not required to analyse the hyperbolicity, the added-mass exchange rate will be approximated by a linear relaxation model

(2.5)\begin{equation} S_a = \frac{\rho_f \alpha_f \alpha_p}{\tau_a} ( c_m^{\star} - c_m ), \end{equation}

with time scale

(2.6)\begin{equation} \tau_a = \frac{4 d_p^2 \alpha_f^{\star}}{3 \nu_f C_D Re_p \alpha_f} . \end{equation}

Physically, $\tau _a$ is the time scale characterizing the expansion/contraction/formation of particle wakes. For example, when a particle moves from a region with large $\alpha _p$ to one with small $\alpha _p$ (i.e. to larger spacing between particles), $c_m$ will be smaller than $c_m^{\star }$. Thus, the wake will grow by entraining fluid with velocity $\boldsymbol {u}_f$ and kinetic energy $u_f^2/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 Reference Moore and Balachandar2019).

2.3. Added-mass function

In principle, by formulating a physically accurate function for $c_m^{\star }$, the two-fluid model will be able to account correctly for unsteady effects. For example, $c_m^{\star }$ might depend on $Ac$ (Odar & Hamilton Reference Odar and Hamilton1964), 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^{\star } (x)$ must have the properties $0 \le \alpha _p c_m^{\star } \le 1$ and $c_m^{\star }(0)=C_m$ where $C_m$ is the added-mass constant, which is equal to 1/2 for a spherical particle when $\alpha _p=0$ (Zuber Reference Zuber1964). In addition, one might require $c_m^{\star } (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. (Reference Sangani, Zhang and Prosperetti1991); however, these expressions are valid for $\alpha _p < 0.5$. From the hyperbolicity analysis in § 3, we find that $0.085 < c_m^{\star } < 1/ \alpha _p$, which corresponds physically to $0 < \alpha _f^{\star } < \alpha _f$. These observations suggest the use of the expression proposed by Zuber (Reference Zuber1964) (written to account for the difference between $\boldsymbol {u}_f$ and $\boldsymbol {v}_f$) of the form

(2.7)\begin{equation} c_m^{\star} = \tfrac{1}{2} \min \left( 1 + 2 \alpha_p , 2 \right) . \end{equation}

Sangani et al. (Reference Sangani, Zhang and Prosperetti1991) 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 default expression in the proposed two-fluid model. Nonetheless, as done in Moore & Balachandar (Reference Moore and Balachandar2019) 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 $\alpha _p^{\star }$ versus $\alpha _p$) is the linear form

(2.8)\begin{equation} c_m^{\star} = C_m + (1-C_m) x, \end{equation}

with $x=\alpha _p$ or $x=\alpha _p^{\star }$, and $0 \le C_m \le 2$. Based on their experimental results, Odar & Hamilton (Reference Odar and Hamilton1964) found that $C_m$ depends on the acceleration number as

(2.9)\begin{equation} C_m = 1 - \tfrac{1}{2} e^{- \beta Ac}, \end{equation}

where $\beta \approx 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 Reference Auton, Hunt and Prud'homme1988; Sangani et al. Reference Sangani, Zhang and Prosperetti1991; Mei & Adrian Reference Mei and Adrian1992) suggest that the decomposition of the virtual-mass and history forces used by Odar & Hamilton (Reference Odar and Hamilton1964) is not reliable, and that $C_m=1/2$ independent of $Ac$. In any case, using (2.8) and given that $\alpha _p^{\star } = \alpha _p + c_m \alpha _p (1- \alpha _p)$, at steady state where $c_m=c_m^{\star }$ with $x=\alpha _p$, the value of $\alpha _p$ is the root in the interval $[0 \ 1]$ of a cubic polynomial for $0 \le C_m \le 2$

(2.10)\begin{equation} (1-C_m) \alpha_p^3 + (2 C_m-1) \alpha_p^2 - (C_m+1) \alpha_p + \alpha_p^{\star} = 0. \end{equation}

For $C_m = 1$, the desired root is $\alpha _f^2 = \alpha _f^{\star }$ (which also holds for (2.7) when $\alpha _f < 1/2$). For other values of $C_m$, the root can be found numerically as illustrated in figure 2.

Figure 2. Steady-state relation between $\alpha _p$ and $\alpha _p^{\star }$ for the added-mass function $c_m^{\star }=C_m+(1-C_m)x$ with three values of $C_m$. $(a)$$x=\alpha _p$. $(b)$$x=\alpha _p^{\star }$. The diagonal line corresponds to $c_m^{\star }=0$. For the function in (2.7), the dependence will be the same as $C_m=1$ for $1/2 < \alpha _p$. Note that the curve for $C_m=1$ is the same for both choices of $x$.

As previously noted, the choice of $c_m^{\star }$ 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 $\alpha _p$, etc., derived from PRDNS, experiments and theory.

2.4. Particle-phase pressure tensor

In fluid–particle flows, the particles have uncorrelated motion due to fluid-mediated interactions and direct collisions (Lhuillier et al. Reference Lhuillier, Chang and Theofanous2013). In recent PRDNS studies of bubbly flow (du Cluzeau, Bois & Toutant Reference du Cluzeau, Bois and Toutant2019; du Cluzeau et al. Reference du Cluzeau, Bois, Toutant and Martinez2020), 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 Reference Biesheuvel and van Wijngaarden1984). Indeed, in kinetic theory the magnitude of the dispersed-phase Reynolds stresses is proportional to the granular temperature $\varTheta _p$. In du Cluzeau et al. (Reference du Cluzeau, Bois, Toutant and Martinez2020), it is found that these terms (referred to as $\boldsymbol {M}^{extra}$ and $\boldsymbol {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 $\partial _{\boldsymbol {x}} \alpha _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 (Reference Batchelor1988) 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)

(2.11)\begin{equation} \partial_x p_p = \partial_x \alpha_p \rho_p \varTheta_p - \frac{\rho_p D}{m_p B} \partial_x \alpha_f, \end{equation}

where $m_p$ is the particle mass. He then used physical reasoning to argue that

(2.12)\begin{equation} \frac{\rho_p D}{m_p B} \propto \rho_f u_{fp}^2 . \end{equation}

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 Reference Fox2019), a dispersion term is also found in $\boldsymbol {F}_{pf}$, but written in terms of $\partial _{\boldsymbol {x}} \rho _f$ and not $\partial _{\boldsymbol {x}} \alpha _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 (Reference Syamlal2011) derived a ‘buoyant-force’ term that extends the fluid pressure in the Archimedes force to include a relative-velocity contribution of the form $\rho _f \alpha _f \boldsymbol {u}_{fp} \otimes \boldsymbol {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 $\boldsymbol {R}$ (which is exact from kinetic theory (Fox Reference Fox2019)) is not present, and the particle-phase pressure tensor term $\partial _{\boldsymbol {x}} \boldsymbol {\cdot } ( \alpha _a \boldsymbol {P}^a_{fp} )$ is replaced by the buoyant-force term $\alpha _p \partial _{\boldsymbol {x}} \boldsymbol {\cdot } ( \rho _f \alpha _f \boldsymbol {u}_{fp} \otimes \boldsymbol {u}_{fp} )$. In § 3.2, we find that the $\boldsymbol {R}$ contribution to the fluid-phase momentum flux is not required for hyperbolicity (and, for constant $\rho _f$, can be combined with the fluid pressure as done in § 5.2). Thus, by rewriting the buoyant-force term in the form

(2.13)\begin{equation} \alpha_p \partial_{\boldsymbol{x}} \boldsymbol{\cdot} ( \rho_f \alpha_f \boldsymbol{u}_{fp} \otimes \boldsymbol{u}_{fp} ) = \partial_{\boldsymbol{x}} \boldsymbol{\cdot} ( \rho_f \alpha_f \alpha_p \boldsymbol{u}_{fp} \otimes \boldsymbol{u}_{fp} ) - \rho_f \alpha_f (\boldsymbol{u}_{fp} \otimes \boldsymbol{u}_{fp} ) \boldsymbol{\cdot} \partial_{\boldsymbol{x}} \alpha_p , \end{equation}

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 (Reference Zhang, Ma and Rauenzahn2006) and Zhang (Reference Zhang2020) 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 $\partial _{\boldsymbol {x}} \boldsymbol {\cdot } ( \alpha _p \boldsymbol {\varSigma }_{pfp} )$, where $\boldsymbol {\varSigma }_{pfp}$ is the PFP stress, appears in their formulation. For uniform potential flow with constant $\rho _f$, Zhang (Reference Zhang2020) finds that

(2.14)\begin{equation} \alpha_p \boldsymbol{\varSigma}_{pfp} = \rho_f \left[ C_1 ( \alpha_p ) u_{fp}^2 \boldsymbol{I} + C_2 ( \alpha_p ) \boldsymbol{u}_{fp} \otimes \boldsymbol{u}_{fp} \right], \end{equation}

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 $\boldsymbol {\varSigma }_{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 $\boldsymbol {\varSigma }_{pfp} \propto \rho _f u_{fp}^2$).

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

(2.15)\begin{equation} \boldsymbol{P}_p = p_p \boldsymbol{I} + \alpha_a \boldsymbol{P}_{fp}^a, \end{equation}

where $p_p = \alpha _p p_p^k + \alpha _a p_f^a$ with $p_p^k = \rho _p \varTheta _p ( 1 + 4 \alpha _p^{\star } g_{0} )$, $p_f^a = \rho _f \varTheta _p ( 1 + 4 \alpha _p^{\star } g_{0} )$ and

(2.16)\begin{equation} \boldsymbol{P}_{fp}^a = \rho_f \frac{2}{3 \gamma_p} \left( \frac{1}{2} u_{fp}^2 \boldsymbol{I} + \boldsymbol{u}_{fp} \otimes \boldsymbol{u}_{fp} \right) \left( 1 + 4 \alpha_p^{\star} g_{0} \right) , \end{equation}

which is a particular form of (2.14). This model for $\boldsymbol {P}_p$ combines the kinetic-theory dependence on $\varTheta _p$ due to uncorrelated velocity fluctuations and direct collisions when $\rho _f \ll \rho _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 Reference van Wijngaarden1976; Batchelor Reference Batchelor1988; Zhang et al. Reference Zhang, Ma and Rauenzahn2006; Zhang Reference Zhang2020, for detailed discussions).

In (2.16), the contribution $1 + 4 \alpha _p^{\star } g_{0}$ with $1 \le 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 $\rho _f \approx \rho _p$ or for deformable particles such as bubbles). For example, one might consider replacing $\alpha _p^{\star } g_{0}$ with $\alpha _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 $\gamma _p$ is not so large as to make $\boldsymbol {P}_{fp}^a$ negligible. Furthermore, in § 3 we find that with $\alpha _p=0$, the system is hyperbolic even when $c_m=0$ (i.e. $\alpha _a=0$ in (2.15)). Thus, it is possible for $\boldsymbol {P}_{fp}^a$ in (2.16) to depend linearly on $\alpha _p^{\star }$ (so that the fluid-mediated pressure depends on $\alpha _p^2$) 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(\alpha _a \boldsymbol {P}_{fp}^a) \propto (\alpha _p^{\star })^2 u_{fp}^2$ 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 (Reference Fox2019), it was shown that the mixture pressure tensor has the form $\boldsymbol {P}_{mix} = \boldsymbol {P}_1 + \boldsymbol {P}_2 + c_{12} \boldsymbol {P}_{BE}$ regardless of the size ratio between the hard spheres. The Boltzmann–Enskog contribution $\boldsymbol {P}_{BE}$ leads to the term involving $\boldsymbol {R}$ in the fluid-phase momentum balance when added mass is neglected ($\alpha _a=0$). Biesheuvel & van Wijngaarden (Reference Biesheuvel and van Wijngaarden1984) derive a contribution to the mixture stress and liquid-phase Reynolds stresses with the same tensorial form as $\boldsymbol {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 $\alpha _a \boldsymbol {P}_{fp}^a$ 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 Reference Fox2019). Finally, note that the contribution $\partial _{\boldsymbol {x}} \boldsymbol {\cdot } ( \alpha _a \boldsymbol {P}_{fp}^a )$ 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 $\gamma _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(\boldsymbol {P}_{fp}^a) \propto 5/(3 \gamma _p)$). On the other hand, the tensorial form of $\boldsymbol {P}_{fp}^a$ must be kept consistent with that of $\boldsymbol {R}$ as both arise due from the same term in the kinetic-theory derivation (Fox Reference Fox2019). As seen from (2.14), up to the scalar coefficients that can depend on $\alpha _p$, this tensorial form is the only one that can be formed from $\boldsymbol {u}_{fp}$ (Zhang Reference Zhang2020).

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 $\rho _p \ll \rho _f$. This is consistent with Batchelor (Reference Batchelor1988) 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. (Reference Zhang, Ma and Rauenzahn2006); Zhang (Reference Zhang2020) who demonstrate that an inter-species stress arises due to particle–fluid–particle interactions, which do not depend on $\varTheta _p$. Nonetheless, future research should be devoted to refining the model for $\alpha _a \boldsymbol {P}_{fp}^a$ in (2.16) to account for the $\alpha _p$ and $Re_p$ dependencies observed in PRDNS.

2.5. Limiting cases

Having previously investigated the case with $\rho _f \ll \rho _p$ (Fox Reference Fox2019), in the remainder of this work we are particularly interested in the limiting cases with $\rho _p=0$ (i.e. the particles have zero mass) so that $\rho _e \alpha _p^{\star } = \rho _f \alpha _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.

3.1. 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

(3.1)\begin{equation} \boldsymbol{X} := \left( X_1 , X_2, X_3, X_4, X_5, X_6, X_7, X_8 \right)^t = \left(\alpha_p , \rho_e \alpha_p^{\star} / \rho_p , Z \alpha_f^{\star} , u_p, u_f, E_p, E_f, k_f \right)^t . \end{equation}

We define $Z$ such that $\rho _f = Z \rho _p$. As $\rho _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 $( \alpha _p, Z, \alpha _f^{\star } )$ by

(3.2a–c)\begin{equation} \alpha_p = X_1 , \quad Z = \frac{X_3 + X_2 - X_1}{\alpha_f}, \quad \alpha_f^{\star} = \frac{X_3}{Z} , \end{equation}

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. $\boldsymbol {R}$ and $\boldsymbol {r}$ in the fluid phase) and forces (i.e. $\boldsymbol {F}_{fp}$ and $D_{fp}$) that are specific to hard-sphere mixtures (Fox Reference Fox2019).

Table 3. One-dimensional compressible two-fluid model with the densities, pressures, fluxes and forces normalized by $\rho _p$. The reference pressure $p^{\star }_o$ is constant and has the same units as $\varTheta _f$, and $Z_0$ is the reference density ratio.

Table 4. Simplified version of 1-D compressible two-fluid model from table 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.

Formally, the 1-D models can be rewritten as

(3.3)\begin{equation} \boldsymbol{A} ( \boldsymbol{X}) \partial_t \boldsymbol{X} + \boldsymbol{B} (\boldsymbol{X}) \partial_x \boldsymbol{X} = \boldsymbol{0}, \end{equation}

where the coefficient matrices $\boldsymbol {A}$ and $\boldsymbol {B} := \boldsymbol {B}^* + \boldsymbol {B}_0$ yield the flux matrix $\boldsymbol {F} := \boldsymbol {A}^{-1} \boldsymbol {B}$. Here, $\boldsymbol {B}_0$ is the contribution due to the pressure and buoyancy terms, and $\boldsymbol {A}$ and $\boldsymbol {B}^*$ are from the convection terms. Written in terms of the components of $\boldsymbol {X}$, the latter are

(3.4)\begin{equation} \boldsymbol{A} := \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & X_2 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & X_3 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & X_2 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & X_3 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix} \end{equation}

and

(3.5)\begin{equation} \boldsymbol{B}^* := \begin{bmatrix} X_4 & 0 & 0 & X_1 & 0 & 0 & 0 & 0 \\ 0 & X_4 & 0 & X_2 & 0 & 0 & 0 & 0 \\ 0 & 0 & X_5 & 0 & X_3 & 0 & 0 & 0 \\ 0 & 0 & 0 & X_2 X_4 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & X_3 X_5 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & X_2 X_4 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & X_3 X_5 & 0 \\ 0 & 0 & 0 & 0 & 2X_8 & 0 & 0 & X_5 \end{bmatrix} . \end{equation}

The components of $\boldsymbol {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 $\boldsymbol {F}$ can be written $u_f + u_0 \lambda _k$ with $k \in \{1,\ldots ,8\}$, where for fixed values of $\overline {p_f^o} = {Z_0 p^{\star }_o}/{u_0^2} = {p_f^o}/{\rho _p u_0^2}$, $\gamma _f$ and $\gamma _p$, each $\lambda _k$, called here a normalized eigenvalue, depends on five dimensionless parameters

(3.6a–e)\begin{equation} \alpha_p, \quad c_m, \quad Ma_s = \frac{u_p}{u_{0}}, \quad \varTheta_r = \frac{\varTheta_p}{u_{0}^2}, \quad K_r = \frac{k_f}{u_{0}^2}, \end{equation}

where $\varTheta _f = u_{0}^2 + \varTheta _{0}$ and $\varTheta _{0}$ is defined by $p_f=0$ from the stiffened-gas model. The parameter $\varTheta _0$ depends on $Z$. The $\lambda _k$ are the roots of $P(X)=Q(u_0X)/u_0^8$, where $Q$ is the characteristic polynomial of $\boldsymbol {F}-u_f \boldsymbol {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 = (\gamma _f p_f^o )^{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 $\alpha _p=0$, the normalized eigenvalues (which are the same for both models) can be computed analytically

(3.7a–d)\begin{align} & 0, \quad \pm \sqrt{\gamma_f + \gamma_f ( \gamma_f - 1) \frac{\overline{p_f^o}}{Z} + 6 K_r} , \quad Ma_s,\nonumber\\ &\quad \frac{1+ (1+ 1/\gamma_p)c_m Z }{1 + c_m Z} Ma_s \pm \sqrt{ \frac{( 1 + (1+ 1/ \gamma_p^2) c_m Z ) c_m Z }{(1+c_m Z)^2} Ma_s^2 + \gamma_p \varTheta_r } \end{align}

and are always real valued, including when $c_m=0$. Here, the two ‘fluid-phase’ eigenvalues that depend on $\overline {p_f^o}$ are always real and distinct. Note that when $\varTheta _r=0$ the ‘particle-phase’ eigenvalues scale with $Ma_s$, but always remain real-valued. When $Ma_s=0$, these eigenvalues depend on $\varTheta _r$ like an ideal gas ($\gamma _p=5/3$). In the following, we investigate the behaviour of the eigenvalues for a stiffened gas ($\gamma _f=29/4$) with fixed values of $Z$, namely, $+\infty$, 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 $\alpha _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 $\alpha _p$ mainly due to $g_0$, while their values at $\alpha _p =0$ depend on $c_m$ as shown in (3.7a–d).

Figure 3. Normalized eigenvalues dependent on $\alpha _p$ for the full $(a{,}c{,}e)$ and simplified $(b{,}d{,}f)$ 1-D models. The two eigenvalues dependant on $p_o^{\star }=10^{8}$, and the eigenvalue at $Ma_s$ for the simplified model, are not shown. All eigenvalues are real valued with these parameters as shown in § 3.2. $(a,\!b)\,\, Z = 0.0001, c_{m} = 0.5; (c,\!d)\,\, Z = 1,$$c_{m} = 0.5; (e,f)\,\, Z = 10\,000, c_{m} = 0.5.$

3.2. 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. (Reference Chalons, Fox, Laurent, Massot and Vié2017), Sturm's theorem (Sturm Reference Sturm1829) 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 \in \{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 \in \{0,1,\ldots ,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 $\overline {p_f^o}$. Thus, a Taylor expansion of the $S_n$ can be done when $\epsilon =1/\overline {p_f^o}$ tends to zero, for $\gamma _f=29/4$ and $\gamma _p=5/3$. Then, $S_0=1$, $S_1=6$ and the limit of $\epsilon S_n$ when $\epsilon$ tends to zero is studied

(3.8)\begin{gather} \epsilon S_2 \rightarrow \frac{145}{8} \frac{\alpha_f(1+c_mZ) + \alpha_pZ }{\alpha_fZ(\alpha_f c_m Z + 1)( 1 - c_m\alpha_p)}, \end{gather}

(3.9)\begin{gather} \epsilon S_3 \rightarrow \frac{2175}{16} \frac{\alpha_f + Z\alpha_p + \alpha_f c_m Z}{\alpha_f Z(\alpha_f c_m Z + 1)(1 - c_m\alpha_p)}, \end{gather}

(3.10)\begin{align} & \epsilon S_4 \rightarrow \varTheta_r\frac{p_2(\alpha_f,c_m,Z)}{1 - c_m\alpha_p} +Ma_s^2\frac{p_1(\alpha_f,c_m,Z)}{(1 - c_m\alpha_p)^3} \left\{5 \alpha_f^8(1 - c_m\alpha_p)^2 + c_m Z^2 q_1(\alpha_f,c_m) \right.\nonumber\\ &\quad \left. + Z(1 - c_m\alpha_p)[18\alpha_f^5\alpha_pc_m^2(1+\alpha_f)(1 - c_m\alpha_p)+q_2(\alpha_f,c_m)] \right\}, \end{align}

(3.11)\begin{align} & \epsilon S_5 \rightarrow \frac{\mu(\alpha_p,c_m,Z,Ma_s, \varTheta_r)^2 } {(1 - c_m\alpha_p)} \left\{ \varTheta_r p_3(\alpha_f,c_m,Z) +Ma_s^2\left[ c_m Z^2 q_4(\alpha_f,c_m) \right.\right.\nonumber\\ &\left.\left. \quad +5(1+\alpha_fc_mZ)[q_3(\alpha_f,c_m)+6c_m^2\alpha_f^5\alpha_p(1-\alpha_f)(1 - c_m\alpha_p)] \right] \right\}, \end{align}

where each $p_k(\alpha _f,c_m,Z)$ is a polynomial function of $\alpha _f$, $c_m$ and $Z$ that is positive for $\alpha _f\in [0,1]$, $c_m\ge 0$ and $Z\ge 0$, each $q_k(\alpha _f,c_m)$ is a polynomial function of $\alpha _f$ and $c_m$ that is positive for $\alpha _f\in [0,1]$ if $c_m$ is large enough, a sufficient condition being $c_m>0.085$. Moreover, $\mu$ is a function of the parameters $(\alpha _p,c_m,Z,Ma_s, \varTheta _r)$. Then the limits of the $\epsilon S_k$ are positive as long as $1 - c_m\alpha _p>0$, and if neither $c_m$ nor $\varTheta _r$ are too small, $c_m > 0.085$ being a sufficient condition. The simplified model is then hyperbolic in those cases in the limit of infinite $\overline {p_f^o}$.

3.3. Eigenvalues for specific cases

We are specifically interested to know whether the 1-D models are hyperbolic for all physically relevant values of $0 \le \alpha _p \le 1$ and $0 \le c_m \le 1/\alpha _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, $\varTheta _r$ mainly affects the ‘particle-phase’ eigenvalues and $\varTheta _r=0$ is known to yield complex eigenvalues when added mass is neglected (Fox Reference Fox2019); hence, the analysis of this limiting case is of particular interest. As done in Fox (Reference Fox2019), we will make use of stability plots found from the Sturm coefficients to check for complex eigenvalues in $\alpha _p$–$c_m$ parameter space.

3.3.1. Granular flow with $\rho _f=0$

The limit $Z=0$ corresponds to a granular flow with $\rho _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 $\varTheta _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 $\varTheta _r=0$.

Figure 4. Hyperbolicity plot for the full $(a{,}b{,}e{,}f)$ and simplified $(c{,}d{,}g{,}h)$ 1-D models for granular flow ($Z = 0$) with varying $Ma_s$. The Sturm test function is negative in black regions, which correspond to unphysical values of $c_m$ as discussed in § 3.2. The minimum value of $\varTheta _r$ needed to avoid round-off error is shown. $(a)\, Ma_{s} = 0, \varTheta_{r} = 0;$$(b)\, Ma_{s} = 0.1, \varTheta_{r} = 10^{-11};$$(c)\, Ma_{s} = 0, \varTheta_{r} = 0;$$(d)\, Ma_{s} = 0.1, \varTheta_{r} = 10^{-15};$$(e) \, Ma_{s} = 1, \ \varTheta_{r} = 10^{-9};$$(\,f) \ Ma_{s} = 10, \varTheta_{r} = 10^{-7};$$(g)\, Ma_{s} = 1, \varTheta_{r} = 10^{-14};$$(h)\, Ma_{s} = 10, \varTheta_{r} = 10^{-12}.$

3.3.2. Neutrally buoyant flow with $\rho _p=\rho _f$

The limit $Z=1$ corresponds to a neutrally buoyant flow with $\rho _p=\rho _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 $\alpha _p=0$ can be less that unity for values of $\alpha _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.

Figure 5. Hyperbolicity plot for the full $(a{,}b{,}e{,}f)$ and simplified $(c{,}d{,}g{,}h)$ 1-D models for neutrally buoyant flow ($Z=1$) and varying $Ma_s$. The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues. As shown in the analysis of the Sturm coefficients in § 3.2, only the black regions with $c_m < 0.085$ are of interest. $(a{,}c)\, Ma_{s} = 0; (b{,}d)\, Ma_{s} = 0.1; (e{,}g)\,\, Ma_{s} = 1; (\,f{,}h)\, Ma_{s} = 10.$

Figure 6. Eigenvalues dependent on $\alpha _p$ for the full $(a{,}c{,}e{,}g)$ and simplified $(b{,}d{,}f{,}h)$ 1-D models. Complex eigenvalues are observed in (e–h) (a–d, $c_m=0.08$; e–h, $c_m=0.008$) and only the real parts are plotted. Both models yield similar eigenvalues. The two eigenvalues that become complex correspond to the particle-phase eigenvalues at $\alpha _p=0$ in (3.7a–d). $(a{,}b)\ Z = 1,\ c_m = 0.08;$$(c{,}d)\, Z = 10\,000, c_m = 0.08;$$(e{,}f)\, Z = 1, c_m = 0.008;$$(g{,}h)\, Z = 10\,000,\ c_m = 0.008.$

3.3.3. Bubbly flow with $\rho _p=0$

The limit $Z \to \infty$ corresponds to a bubbly flow with $\rho _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^{\star }$ and $\tau _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.

Figure 7. Hyperbolicity plot for the full $(a{,}b{,}e{,}f)$ and simplified $(c{,}d{,}g{,}h)$ 1-D models for bubbly flow ($Z \to +\infty$) and varying $Ma_s$. The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues. As shown in the analysis of the Sturm coefficients in § 3.2, only the black regions with $c_m < 0.085$ are of interest. $(a{,}c)\, Ma_s = 0; (b{,}d)\, Ma_s = 0.1;\ (e{,}g)\, Ma_s = 1; \ (\,f{,}h)\, Ma_s = 10.$

In conclusion, it is worth noting that in the full model the region of hyperbolicity for $0.7 < \alpha _p$ includes $c_m \to 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 $\alpha _p \approx 0.1$, regardless of the material-density ratio.

4.1. Conservative form of 1-D model

The added-mass contribution to the particle-phase pressure tensor $\boldsymbol {P}_{fp}^a$ is purely mechanical. This implies that the granular energy balance for $\varTheta _p$ has a compression source term that depends only on the ‘thermodynamic’ pressure $p_p$ (see appendix A in Houim & Oran Reference Houim and Oran2016, 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 $\gamma _p = 5/3$)

(4.1)\begin{align} &\tfrac{3}{2} \left( Y_2 \partial_t \varTheta_p + Y_2 X_4 \partial_x \varTheta_p \right) = - p_p \partial_x X_4 - 2 \alpha_p^{\star} \varTheta_p [ (X_4-X_5) \partial_x Z - Z \partial_x X_5 ] \nonumber\\ &\quad + K [ 2 (1-a) X_8 - 3a \varTheta_p ] + \tfrac{1}{2} \max(S_a , 0) \left[ (X_4-X_5)^2 - 3 \varTheta_p \right] \end{align}

which has the necessary property that the right-hand side is non-negative when $\varTheta _p$ is null. In contrast, if $p_p$ were replaced by $(p_p+P_{fp}^a)$ in (4.1), then the first term on the right-hand side would be negative when $\varTheta _p=0$ and $\partial _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 mixture mass ($\varrho =Y_2+Y_3$), momentum $\mathcal {M}=Y_4+Y_5$ and energy $\mathcal {E}=Y_6+Y_7$ balances can be written as

(4.2)\begin{equation} \left. \begin{array}{c@{}} \partial_t \varrho + \partial_x \mathcal{M} = 0 \\ \partial_t \mathcal{M} + \partial_x( Y_4 X_4 + Y_5 X_5 + P_f + p_p + \alpha_p^{\star} Z R ) = \varrho g_x \\ \partial_t \mathcal{E} + \partial_x ( Y_6 X_4 + Y_7 X_5 + \alpha_f^{\star} P_f X_5 + \alpha_p^{\star} P_f X_4 + p_p X_4 + \alpha_p^{\star} Z R X_4 + r ) = \mathcal{M} g_x \end{array} \right\} \end{equation}

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

(4.3)\begin{equation} \left. \begin{array}{c@{}} \partial_t Y_1 + \partial_x ( Y_1 X_4 ) = 0 \\ \partial_t Y_2 + \partial_x Y_4 = S_a \\ \partial_t Y_4 + \partial_x ( Y_4 X_4 + p_p ) = Y_2 g_x - \partial_x ( \alpha_a P_{fp}^a )- \alpha_p^{\star} \partial_x P_f - F_{pf} + K (X_5-X_4) - S_{pf} \\ \partial_t Y_6 + \partial_x ( Y_6 X_4 + p_p X_4 ) = Y_4 g_x - X_4 \partial_x ( \alpha_a P_{fp}^a ) - X_4 \alpha_p^{\star} \partial_x P_f - D_{pf} - D_E - S_E \\ \partial_t Y_8 + \partial_x ( Y_8 X_5 ) = - 2 X_3 X_8 \partial_x X_5 + D_{PT} - X_3 \varepsilon_{PT} \end{array} \right\} \end{equation}

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 $\varTheta _p$ shown in (4.1), which is preferable for dissipative systems to avoid round-off errors due to the very small value of $\varTheta _p$ (see Houim & Oran (Reference Houim and Oran2016), for a discussion of this point), and for the numerical treatment of body forces.

4.2. Numerical solver

The 1-D model in table 5 has the form

(4.4)\begin{equation} \partial_t \boldsymbol{Y} + \partial_x\, \boldsymbol{f}(\boldsymbol{Y}) = \boldsymbol{h} (\boldsymbol{Y}), \end{equation}

where $\boldsymbol {f}$ are the spatial fluxes, and $\boldsymbol {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

(4.5)\begin{equation} \boldsymbol{Y}^{n+1}_i = \boldsymbol{Y}^n_i - \frac{\Delta t}{\Delta x} (\, \boldsymbol{f}^n_{i}-\boldsymbol{f}^n_{i-1} ) + \Delta t \, \boldsymbol{h} (\boldsymbol{Y}^n_{i+1}, \boldsymbol{Y}^n_i, \boldsymbol{Y}^n_{i-1}) . \end{equation}

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_4^0 = Y_6^0=0$ yields $Y_4^1 = \Delta t Y_2^0 g_x$, $Y_6^1=0$). In general, if there is mean slip between the phases (i.e. $Z_0 \neq 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.

The numerical fluxes in (4.5) are defined using a classical HLL approach (Harten, Lax & van Leer Reference Harten, Lax and van Leer1983; Toro Reference Toro1997)

(4.6)\begin{equation} \boldsymbol{f}^n_{i} = \frac{ a^+ \boldsymbol{f} ( \boldsymbol{Y}^n_i ) - a^- \boldsymbol{f}( \boldsymbol{Y}^n_{i+1} ) + a^+ a^- ( \boldsymbol{Y}^n_{i+1} - \boldsymbol{Y}^n_{i} ) }{a^+ - a^-}, \end{equation}

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^- \approx - a^+$. In our simulations, the eigenvalues are computed at each time step from the characteristic polynomial. On a uniform grid with spacing $\Delta x$, the time step is set using $\Delta x / \Delta 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 Reference Toro1997), will generate significant numerical diffusion for the system in table 5. For example, the effective diffusivity of the volume fraction is $D \propto \Delta x \, a^+$ due to the final difference term on the right-hand side of (4.6). For this reason, unless $\Delta 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 Reference Toro1997).

4.3. 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 $\nu _f =10^{-6}\ \textrm {m}^{2}\,\textrm {s}^{-1}$. The fluid temperature $\varTheta _f$ is set to be $5000\ \textrm {m}^{2}\,\textrm {s}^{-2}$ larger than $\varTheta _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}\ \textrm {m}$ and set $c_m^{\star } = 1/2$. Finally, in order to keep the time step reasonable, we set $p_o^{\star } = 10^4\ \textrm {m}^{2}\,\textrm{s}^{-2}$, which yields $a^+ \approx 775\ \textrm {m}\,\textrm {s}^{-1}$. Increasing $p_o^{\star }$ will result in smaller variations in $Z$, but requires a correspondingly smaller time step. The computational domain is taken as $x \in [-1/2,1/2]\ \textrm {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 $\Delta x = 1/N\ \textrm {m}$ with $N=(1000,\ 2000,\ 4000)$.

In the first example, we consider a case with $Z_0=1$. The initial conditions are $\alpha _p=0$ on the left half and $\alpha _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 $\alpha _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 $\alpha _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 $\varTheta _p$. For the exact solution, $\varTheta _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.

Figure 8. Primitive variables at $t=0.1\ \textrm {s}$ for Riemann problem with $Z_0=1$ and $\Delta x = 1/N\ \textrm {m}$ (blue, $N=1000$; red, $N=2000$; gold, $N=4000$). The exact solution for $\alpha _p$ is a step function at $x=0$. Here, the HLL fluxes result in numerical diffusion of the volume fractions, which can be reduced by increasing $N$.

In the second example, we consider a case with $Z_0=10^4$ corresponding to buoyant particles. The initial conditions are again $\alpha _p=0$ on the left half and $\alpha _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, $\alpha _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.

Figure 9. Primitive variables at $t=0.1\ \textrm {s}$ with $Z_0=10^4$ and $\Delta x = 1/N\ \textrm {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 $\alpha _p$.

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, $\varTheta _p$ (and, hence, $p_p$) will be much smaller due to inelastic collisions and lubrication effects (Abbas et al. Reference Abbas, Climent, Parmentier and Simonin2010). 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.