2.1 Boltzmann–Enskog kinetic equations

We begin by introducing the VDFs $f_{\unicode[STIX]{x1D6FC}}(t,\boldsymbol{x},\boldsymbol{v})$ for each type of particle $\unicode[STIX]{x1D6FC}=1,2$ with spatial location $\boldsymbol{x}\in \mathbb{R}^{3}$ , and velocity $\boldsymbol{v}\in \mathbb{R}^{3}$ . The integration of the VDFs over velocity phase space yields the total mass of particles of type $\unicode[STIX]{x1D6FC}$ per unit volume (i.e. the mass concentration), which we denote by $\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}(t,\boldsymbol{x})$ . These VDFs are assumed to obey the following coupled BE kinetic equations (Chapman & Cowling Reference Chapman and Cowling1952; Tham & Gubbins Reference Tham and Gubbins1971; Barajas et al. Reference Barajas, Garcia-Colin and Pinã1973; Marchisio & Fox Reference Marchisio and Fox2013) (the same equations are used with the Boltzmann collision operator (Grad Reference Grad and Devienne1960; Sirovich Reference Sirovich1962; Hamel Reference Hamel1966), but with a simpler definition for ${\mathcal{C}}_{\unicode[STIX]{x1D6FC}}$ . The $N$ -particle probability function used in Zhang et al. (Reference Zhang, Ma and Rauenzahn2006) can be used to derive (2.1) for the one-particle VDF for binary interactions. Nonetheless, the closure of the collision terms requires a model for the two-particle VDF expressed in terms of the one-particle VDF)

(2.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\,f_{1}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{x}}\,f_{1}={\mathcal{C}}_{1},\quad \unicode[STIX]{x2202}_{t}\,f_{2}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{x}}\,f_{2}={\mathcal{C}}_{2}.\end{eqnarray}$$

The first two terms on the left-hand sides are accumulation and free transport, respectively. The terms ${\mathcal{C}}_{\unicode[STIX]{x1D6FC}}$ on the right-hand sides correspond to particle–particle collisions, and include contributions due to like and unlike particles. In the model equations below, all collisions are assumed to be elastic so that the total energy of the system is conserved. The exact forms of the (inelastic) collision terms are given in appendix A, along with the derivation of the collisional contributions to the moment equations starting from (2.1). The inelastic forms of the collision terms can be used to describe granular flows, but a thermal energy balance would then be required to account for the loss of kinetic energy due to collisions. (Alternatively, in the fluid–particle limit, the kinetic energy dissipated in particle–particle collisions could be included as a source term in the fluid energy balance.) Results for the moment equations with the Boltzmann collision operator can be found in the literature (Grad Reference Grad and Devienne1960; Sirovich Reference Sirovich1962; Hamel Reference Hamel1966). Here, it suffices to remark that the BE model accounts for the finite size of the particles and these additional terms are the primary focus of the present work. For clarity, in the moment equations these terms will be identified using the subscript $\text{}_{BE}$ .

2.2 Definition of hard-sphere model parameters

In the following sections, we consider the case where $d_{1}\leqslant d_{2}$ (and eventually $d_{1}\ll d_{2}$ , corresponding to the fluid–particle limit). Hereinafter, $V_{\unicode[STIX]{x1D6FC}}=(\unicode[STIX]{x03C0}/6)d_{\unicode[STIX]{x1D6FC}}^{3}$ and $m_{\unicode[STIX]{x1D6FC}}$ are the volume and mass, respectively, of a single particle of type $\unicode[STIX]{x1D6FC}$ . Note that the total volume of particles of type $\unicode[STIX]{x1D6FC}$ per unit volume is $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FC}}:=(V_{\unicode[STIX]{x1D6FC}}/m_{\unicode[STIX]{x1D6FC}})\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}$ . Thus, $0<\unicode[STIX]{x1D711}_{1}+\unicode[STIX]{x1D711}_{2}\,<\,1$ , i.e. the sum of the particle volume concentrations do not sum to unity. The particles are assumed to have no internal degrees of freedom and, hence, (2.1) will conserve the total mass, momentum and energy of the particles. Body forces are not included in (2.1) (see appendix B for details on how to add them) because they do not change the spatial fluxes. The definitions of the other parameters appearing in the kinetic model are as follows. The collisional cross-section diameter for unlike collisions is

(2.2) $$\begin{eqnarray}d:={\textstyle \frac{1}{2}}(d_{1}+d_{2}),\end{eqnarray}$$

from which we define the geometric parameter $\unicode[STIX]{x1D712}_{\unicode[STIX]{x1D6FC}}:=d/d_{\unicode[STIX]{x1D6FC}}$ . The mass ratios are defined by

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}:=\frac{2m_{\unicode[STIX]{x1D6FC}}}{m_{1}+m_{2}},\end{eqnarray}$$

and the material densities of the particles by $\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}$ . Because these densities are assumed constant, we have $m_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D70C}_{\unicode[STIX]{x1D6FC}}V_{\unicode[STIX]{x1D6FC}}$ . By definition,

(2.4a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{1}+\unicode[STIX]{x1D707}_{2}=2,\quad \text{and}\quad \unicode[STIX]{x1D70C}^{\star }:=\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D712}_{2}^{3}\unicode[STIX]{x1D707}_{2}=\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D712}_{1}^{3}\unicode[STIX]{x1D707}_{1}={\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D70C}_{2}(d_{1}+d_{2})^{3}}{4(\unicode[STIX]{x1D70C}_{1}d_{1}^{3}+\unicode[STIX]{x1D70C}_{2}d_{2}^{3})}},\end{eqnarray}$$

which can be used to demonstrate the conservation of momentum and energy. The mean velocities $\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}$ and velocity covariance tensors $\unicode[STIX]{x1D748}_{\unicode[STIX]{x1D6FC}}$ are found from $f_{\unicode[STIX]{x1D6FC}}$ as described in § 3. The phase temperatures are then defined in terms of the trace $\unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D6FC}}:=\frac{1}{3}\text{tr}(\unicode[STIX]{x1D748}_{\unicode[STIX]{x1D6FC}})$ , and the ‘mixture’ temperature by

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6E9}:=\unicode[STIX]{x1D6E9}_{1}+\unicode[STIX]{x1D6E9}_{2}+{\textstyle \frac{1}{3}}|\boldsymbol{u}_{12}|^{2},\end{eqnarray}$$

where the mean-slip velocity is defined by $\boldsymbol{u}_{12}:=\boldsymbol{u}_{1}-\boldsymbol{u}_{2}=-\boldsymbol{u}_{21}$ .

2.3 Definition of particle-pair distribution function

In the BE kinetic theory of hard-sphere mixtures with finite-size particles, the particle-pair distribution function appearing in the collision integrals (see appendix A for more details) plays a key role as it determines which forces and fluxes beyond the perfect-gas description will appear in the two-fluid model. Thus, in order to produce the Archimedes force in the fluid–particle limit, we propose the following approximation for the particle-pair distribution function with the geometric parameter $\frac{1}{2}\leqslant y\leqslant 1$ (aside from the factors $h_{\unicode[STIX]{x1D6FC}}$ and ${\mathcal{V}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ , this model is essentially the same as the truncated series expansion for dense gases proposed by Enskog (Reference Enskog1911a ); see chap. 16 in Chapman & Cowling (Reference Chapman and Cowling1952) for details)

(2.6) $$\begin{eqnarray}\displaystyle & & \displaystyle f^{(2)}(\boldsymbol{x}-yd\,\boldsymbol{x}_{12},m_{1},\boldsymbol{v}_{1};\boldsymbol{x}+(1-y)d\,\boldsymbol{x}_{12},m_{2},\boldsymbol{v}_{2})\nonumber\\ \displaystyle & & \displaystyle \quad \approx g_{12}\,f_{1}\,f_{2}+{\displaystyle \frac{d}{{\mathcal{V}}_{12}}}\,\boldsymbol{x}_{12}\boldsymbol{\cdot }[(1-y)\unicode[STIX]{x1D71A}_{1}\,f_{1}^{\star }\unicode[STIX]{x2202}_{\boldsymbol{x}}(h_{2}\,f_{2})-y\unicode[STIX]{x1D71A}_{2}\,f_{2}^{\star }\unicode[STIX]{x2202}_{\boldsymbol{x}}(h_{1}\,f_{1})],\end{eqnarray}$$

where the volume ratio is defined by ( ${\mathcal{V}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FC}}=4$ )

(2.7) $$\begin{eqnarray}{\mathcal{V}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}:={\displaystyle \frac{(d_{\unicode[STIX]{x1D6FC}}+d_{\unicode[STIX]{x1D6FD}})^{3}}{d_{\unicode[STIX]{x1D6FC}}^{3}+d_{\unicode[STIX]{x1D6FD}}^{3}}}={\displaystyle \frac{\text{collision volume}}{\text{sum of particle volumes}}}.\end{eqnarray}$$

In (2.6), $f_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x},\boldsymbol{v}_{\unicode[STIX]{x1D6FC}})=\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}(\boldsymbol{x})f_{\unicode[STIX]{x1D6FC}}^{\star }(\boldsymbol{x},\boldsymbol{v}_{\unicode[STIX]{x1D6FC}})$ defines the normalized VDF $f_{\unicode[STIX]{x1D6FC}}^{\star }$ and the unit vector $\boldsymbol{x}_{12}$ between particle centres is defined in (A 1). The model in (2.6) is mathematically consistent in the sense that it treats both types of particles in the same manner, i.e. it obeys the indifferentiability principle (Garzó, Santos & Brey Reference Garzó, Santos and Brey1989; Andries, Aoki & Perthame Reference Andries, Aoki and Perthame2002). When applied to like particles (i.e. $m_{1}=m_{2}$ and $d_{1}=d_{2}$ ), the gradient terms reduce to expressions involving $\boldsymbol{x}_{12}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{x}}h_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FC}}$ . The indifferentiability principle also holds for the collision integrals in appendix A, and is thus inherited from the hard-sphere collision operator.

The first term on the right-hand side of (2.6) is the classical Boltzmann approximation where spatial gradients are neglected. The second term on the right-hand side accounts for spatial gradients at the particle scale. In particular, in the fluid limit where $y\rightarrow 1$ , equation (2.6) models how the spatial gradients in fluid properties vary on the scale of the particle diameter $d_{2}$ . Obviously, equation (2.6) neglects higher-order terms in the expansion about the point $\boldsymbol{x}$ . However, this is consistent with the observations in Barajas et al. (Reference Barajas, Garcia-Colin and Pinã1973). Nonetheless, we must emphasize that the form of the first-order derivative term in (2.6) is chosen to capture the limit where $y\rightarrow 1$ and, hence, it may not be physically accurate for equal-size particles where $y=1/2$ , but with different densities.

As will be discussed in § 3.6, the functions $h_{1}$ and $h_{2}$ are chosen such that the isotropic flux pressure and the buoyancy pressure are identical for each phase. This specification yields

(2.8a,b ) $$\begin{eqnarray}h_{1}:=1+{\mathcal{V}}_{11}\unicode[STIX]{x1D711}_{1}g_{0}+{\mathcal{V}}_{12}\unicode[STIX]{x1D711}_{2}g_{12},\quad h_{2}:=1+{\mathcal{V}}_{22}\unicode[STIX]{x1D711}_{2}g_{0}+{\mathcal{V}}_{21}\unicode[STIX]{x1D711}_{1}g_{12}.\end{eqnarray}$$

In order to arrive at this definition for $h_{\unicode[STIX]{x1D6FC}}$ , the expansion point in the BE collision integral, defined by the choice of $y$ , must be chosen such that the indifferentiability principle holds (these definitions provide alternative expressions for $\unicode[STIX]{x1D70C}^{\star }=({\mathcal{V}}_{12}\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D707}_{2}/8y)=({\mathcal{V}}_{12}\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D707}_{1}/8(1-y))$ , which we use to eliminate $\unicode[STIX]{x1D70C}^{\star }$ and $y$ in the moment equations in § 3. Note that specifying $y$ fixes the interspecies stress tensor contribution in both phases, and thus overcomes the symmetry issue discussed in Zhang et al. (Reference Zhang, Ma and Rauenzahn2006))

(2.9) $$\begin{eqnarray}y:={\displaystyle \frac{d_{2}^{3}}{d_{1}^{3}+d_{2}^{3}}}={\displaystyle \frac{\text{volume of particle 2}}{\text{sum of particle volumes}}}\quad \Longrightarrow \quad {\mathcal{V}}_{12}=8y\unicode[STIX]{x1D712}_{2}^{3}=8(1-y)\unicode[STIX]{x1D712}_{1}^{3}={\mathcal{V}}_{21}.\end{eqnarray}$$

Note that $y$ depends only on the particle diameters and the $h_{\unicode[STIX]{x1D6FC}}$ are symmetric with respect to $\unicode[STIX]{x1D6FC}=1,2$ . Furthermore, ${\mathcal{V}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ , $y$ and $h_{\unicode[STIX]{x1D6FC}}$ are purely geometrical and do not depend on the particle masses. The first contribution to $h_{\unicode[STIX]{x1D6FC}}$ is from the Boltzmann approximation for a perfect gas. The second is the Enskog correction due to the excluded volume from particles of the same type, while the third arises from the excluded volume of the other particles. For equal-size particles, equation (2.8) yields $h_{1}=h_{2}=1+4\unicode[STIX]{x1D711}g_{0}$ (i.e. the indifferentiability principle holds); and for a perfect gas (see § 4) $h_{1}=1/(1-\unicode[STIX]{x1D711}_{2})=1/\unicode[STIX]{x1D6FC}_{f}$ so that $\unicode[STIX]{x1D71A}_{1}h_{1}=\unicode[STIX]{x1D70C}_{f}$ is the gas-phase density (which should not be confused with $\unicode[STIX]{x1D70C}_{1}$ or $\unicode[STIX]{x1D70C}_{2}$ ). For any other choice of $y$ than (2.9), the indifferentiability principle will not be satisfied for equal-sized particles. Thus, equation (2.9) is unique given the assumption of equality between the flux and buoyancy pressures.

In the limit where $d_{1}\ll d_{2}$ , ${\mathcal{V}}_{12}\rightarrow 1$ and $y\rightarrow 1$ , whereas ${\mathcal{V}}_{12}=4$ and $y\rightarrow 1/2$ when $d_{1}=d_{2}$ . In the fluid–particle limit introduced in § 4, (2.6) reduces to

(2.10) $$\begin{eqnarray}f^{(2)}(\boldsymbol{x}-d\,\boldsymbol{x}_{12},m_{1},\boldsymbol{v}_{1};\boldsymbol{x},m_{2},\boldsymbol{v}_{2})\approx g_{12}\,f_{1}\,f_{2}-{\textstyle \frac{1}{2}}d_{2}\,\boldsymbol{x}_{12}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}\,f_{2}^{\star }\unicode[STIX]{x2202}_{\boldsymbol{ x}}(h_{1}\,f_{1}).\end{eqnarray}$$

The first term on the right-hand side of (2.10) is the molecular chaos hypothesis (Cercignani Reference Cercignani1988) with the excluded volume correction, while the second accounts for the spatial gradient of $h_{1}\,f_{1}$ across a particle of diameter $d_{2}$ . (For the Boltzmann collision operator, $f^{(2)}=f_{1}\,f_{2}$ and, hence, $g_{12}$ and the gradient term arise due to the finite particle size.) Note that when evaluating the collision integral $\int _{\mathbb{S}^{+}}$ , the surface $\mathbb{S}^{+}$ corresponds to the surface of type-2 particles in the fluid–particle limit. In other words, the last term in (2.10) models the changes in $h_{1}\,f_{1}$ (i.e. the fluid-particles’ VDF) on the surface of type-2 particles. In principle, additional forces such as virtual mass could be accounted for by modifying (2.6), but this is very challenging to do based on physical arguments (Chliamovitch, Malaspinas & Chopard Reference Chliamovitch, Malaspinas and Chopard2017). Nonetheless, as is well known (Lhuillier et al. Reference Lhuillier, Chang and Theofanous2013), inclusion of the Archimedes force often leads to non-hyperbolic two-fluid models, so we restrict our hyperbolicity analysis to (2.10) in this work.

Furthermore, it is important to note that when applied to like particles of type-2 (i.e. macroscale particles), equation (2.10) completely neglects particle–particle velocity correlations mediated by the fluid phase (sometimes referred to as three-way coupling), which are known to be very important, for example, in dense inertial suspensions (Blanc et al. Reference Blanc, Lemaire, Muenier and Peters2013; Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018). Thus, the resulting two-fluid model will at best be accurate for gas–solid flows where the material density of the solid phase is much larger than that of the fluid phase. In such systems, the forms of the inter- and intra-particle stress tensors are largely determined by particle–particle collisions with the fluid playing a negligible role (e.g. modifying the coefficient of restitution). In principle, the $N$ -particle probability function in Zhang et al. (Reference Zhang, Ma and Rauenzahn2006) overcomes this limitation, but appropriate (closed) forms for the particle–particle interaction force must be postulated to close the moment system arising from the Liouville equation. An alternative approach could be to modify the stress tensors and interaction forces to account for the presence of the fluid phase as suggested in Zhang & Prosperetti (Reference Zhang and Prosperetti1994). Particle-resolved direction-numerical simulations could a useful tool for validating such closures (see, e.g. Tenneti & Subramaniam Reference Tenneti and Subramaniam2014, for examples). However, care must be taken to ensure that the resulting two-fluid model is hyperbolic when developing such closures.

2.4 Definition of radial distribution function

The radial distribution function $g_{12}$ depends on the total volume concentration $\unicode[STIX]{x1D711}:=\unicode[STIX]{x1D711}_{1}+\unicode[STIX]{x1D711}_{2}$ . For a binary, hard-sphere mixture $g_{12}$ can be approximated as (Santos et al. Reference Santos, Yuste and López de Haro1999; Hansen-Goos & Roth Reference Hansen-Goos and Roth2006)

(2.11) $$\begin{eqnarray}g_{12}={\displaystyle \frac{1}{1-\unicode[STIX]{x1D711}}}+\left(g_{0}-{\displaystyle \frac{1}{1-\unicode[STIX]{x1D711}}}\right){\displaystyle \frac{\langle d^{2}\rangle }{\langle d^{3}\rangle }}{\displaystyle \frac{d_{1}d_{2}}{d}},\end{eqnarray}$$

wherein $g_{0}$ is modelled as (Carnahan & Starling Reference Carnahan and Starling1969)

(2.12) $$\begin{eqnarray}g_{0}={\displaystyle \frac{2-\unicode[STIX]{x1D711}}{2(1-\unicode[STIX]{x1D711})^{3}}}.\end{eqnarray}$$

In (2.11), $\langle d^{n}\rangle :=x_{1}d_{1}^{n}+x_{2}d_{2}^{n}$ where $x_{\unicode[STIX]{x1D6FC}}:=(\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}/m_{\unicode[STIX]{x1D6FC}})/(\unicode[STIX]{x1D71A}_{1}/m_{1}+\unicode[STIX]{x1D71A}_{2}/m_{2})$ is the number (or mole) fraction of particles of type $\unicode[STIX]{x1D6FC}$ . This model for $g_{12}$ is a linear interpolation between the limit $d_{1}=d_{2}$ where $g_{12}=g_{0}$ , and the limit $d_{1}\ll d_{2}$ where $g_{12}=(1/1-\unicode[STIX]{x1D711})$ accounts for the excluded volume occupied by the large particles. In Santos et al. (Reference Santos, Yuste and López de Haro1999), the model is validated for diameter ratios down to 0.05 against molecular dynamic simulations. It should be noted that (2.11) does not include non-equilibrium effects due, for example, to a mean-slip velocity in fluid–particle systems, quantified by the mean-slip Mach number $Ma_{s}:=\boldsymbol{u}_{12}/\sqrt{5\unicode[STIX]{x1D6E9}_{1}/3}$ . As discussed in appendix A, it is possible to include corrections to $g_{12}$ that depend on $Ma_{s}$ but further research is needed to validate such corrections.

In this work, we consider only kinetic models with moments up to second order, which is appropriate for collisional flows (Fox Reference Fox2014). The higher-order moments are closed using the anisotropic Gaussian (AG) approximation (Levermore & Morokoff Reference Levermore and Morokoff1996; Vié et al. Reference Vié, Doisneau and Massot2015); most notably the third-order central moments $\unicode[STIX]{x1D64C}$ are null. To go beyond this assumption, it is possible to a derive transport equation for $\unicode[STIX]{x1D64C}$ (and higher-order moments) (Fox et al. Reference Fox, Laurent and Vié2018), or to employ a Chapman–Enskog expansion (Chapman & Cowling Reference Chapman and Cowling1952; Levermore & Morokoff Reference Levermore and Morokoff1996) to express $\unicode[STIX]{x1D64C}$ in terms of the spatial gradients of lower-order moments. However, as such terms do not change the hyperbolicity of the two-fluid model, the AG approximation suffices for our purposes. Finally, we should note that the moment equations derived in § 3 can be applied to polydisperse rapid granular and gas–particle flows with minimal modifications (Kong & Fox Reference Kong and Fox2019).

3.1 Mass balances

The balance equations for mass found from (2.1) are

(3.3a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{1}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}=0,\quad \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{2}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}=0.\end{eqnarray}$$

Here we have used $\int {\mathcal{C}}_{\unicode[STIX]{x1D6FC}}\,\text{d}\boldsymbol{v}=0$ , i.e. the collisions do not change the particle masses. The system density $\unicode[STIX]{x1D71A}:=\unicode[STIX]{x1D71A}_{1}+\unicode[STIX]{x1D71A}_{2}$ is a conserved quantity.

3.2 Momentum balances

The balance equations for mean momentum are

(3.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}\otimes \boldsymbol{u}_{1}+\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D64B}_{1}+\unicode[STIX]{x1D64B}_{12})=\boldsymbol{C}-K\boldsymbol{u}_{12},\\ \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}\otimes \boldsymbol{u}_{2}+\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D748}_{2}+\unicode[STIX]{x1D64B}_{2}+\unicode[STIX]{x1D64B}_{21})=-\boldsymbol{C}+K\boldsymbol{u}_{12}.\end{array}\right\}\end{eqnarray}$$

The collisional pressure tensors (e.g. $\unicode[STIX]{x1D64B}_{1}$ and $\unicode[STIX]{x1D64B}_{12}$ ) and the source terms on the right-hand side of (3.4) are presented below. It is important to recognize that the collisional pressure tensors for like and unlike collisions arise from the same physical mechanisms and, therefore, must both be included in the momentum balances. Moreover, for hard-sphere collisions these pressure tensors are exact (i.e. they do not depend on the form of the velocity VDF). For mono-disperse granular flows, $\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D6FC}}$ is the same as in previous derivations (see, e.g. Jenkins & Savage Reference Jenkins and Savage1983). From (3.4), we can observe that the system momentum $\unicode[STIX]{x1D71A}\boldsymbol{u}:=\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}+\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}$ is a conserved quantity. As discussed in Zhang et al. (Reference Zhang, Ma and Rauenzahn2006), the interspecies stress $\unicode[STIX]{x1D748}^{(12)}=-(\unicode[STIX]{x1D64B}_{12}+\unicode[STIX]{x1D64B}_{21})/\unicode[STIX]{x1D711}_{1}$ is non-zero and contributes to the total mixture stress, as it should to satisfy the indifferentiability principle. However, the expression for $\unicode[STIX]{x1D64F}_{C}$ ( $=\unicode[STIX]{x1D748}^{(21)}$ ) in (10.9) of Zhang & Prosperetti (Reference Zhang and Prosperetti1994) has positive trace, which implies a negative contribution to the mixture pressure. In contrast, the trace of $\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$ is always non-negative, implying a positive contribution to the mixture pressure. In terms of the hyperbolicity of the system, a negative pressure contribution increases the likelihood of having complex eigenvalues.

The drag coefficient on the right-hand side of (3.4), which results directly from the unclosed collision integral in (A 30), is approximated by (Gross & Krook Reference Gross and Krook1956; Morse Reference Morse1963; Andries et al. Reference Andries, Aoki and Perthame2002; Marchisio & Fox Reference Marchisio and Fox2013)

(3.5) $$\begin{eqnarray}K=\unicode[STIX]{x03C0}g_{12}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D71A}_{2}\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D6E9}}{2}}}{\displaystyle \frac{(d_{1}+d_{2})^{2}}{m_{1}+m_{2}}}=6g_{12}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D71A}_{2}\sqrt{{\displaystyle \frac{\unicode[STIX]{x1D6E9}}{2}}}{\displaystyle \frac{(d_{1}+d_{2})^{2}}{\unicode[STIX]{x1D70C}_{1}d_{1}^{3}+\unicode[STIX]{x1D70C}_{2}d_{2}^{3}}}.\end{eqnarray}$$

As drag does not contribute to the spatial fluxes, the exact form of $K$ will not be important when determining whether the two-fluid model is hyperbolic. However, $g_{12}$ and $g_{0}$ appear in the collisional pressure tensors below, so it is important to keep in mind that the model in (2.11) is limited to small $Ma_{s}$ .

In contrast to the drag term, the collisional pressure tensors for like ( $\unicode[STIX]{x1D64B}_{1}$ , $\unicode[STIX]{x1D64B}_{2}$ ) and unlike ( $\unicode[STIX]{x1D64B}_{12}$ , $\unicode[STIX]{x1D64B}_{21}$ ) particles directly affect the system eigenvalues. As shown in appendix A, these symmetric second-order tensors can be computed exactly and are defined by

(3.6) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D64B}_{1}:={\textstyle \frac{1}{5}}{\mathcal{V}}_{11}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D711}_{1}g_{0}(3\unicode[STIX]{x1D6E9}_{1}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D748}_{1}),\quad \unicode[STIX]{x1D64B}_{2}:={\textstyle \frac{1}{5}}{\mathcal{V}}_{22}\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D711}_{2}g_{0}(3\unicode[STIX]{x1D6E9}_{2}\unicode[STIX]{x1D644}+2\unicode[STIX]{x1D748}_{2}),\\ \unicode[STIX]{x1D64B}_{12}:={\textstyle \frac{1}{5}}{\mathcal{V}}_{12}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D711}_{2}g_{12}\unicode[STIX]{x1D64B}_{BE},\quad \unicode[STIX]{x1D64B}_{21}:={\textstyle \frac{1}{5}}{\mathcal{V}}_{12}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D711}_{1}g_{12}\unicode[STIX]{x1D64B}_{BE},\end{array}\right\}\end{eqnarray}$$

where the BE contribution from unlike collisions is

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D64B}_{BE}:={\textstyle \frac{3}{2}}\unicode[STIX]{x1D6E9}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D748}_{2}+\boldsymbol{u}_{12}\otimes \boldsymbol{u}_{12}.\end{eqnarray}$$

As shown in appendix A, the relations in (3.6) and (3.7) are exact in the sense that they do not depend on the shapes of the VDFs. The same is true for the momentum-exchange terms on the right-hand side of (3.4), defined by

(3.8) $$\begin{eqnarray}\boldsymbol{C}:={\textstyle \frac{1}{5}}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D711}_{2}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(2)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}h_{1}\unicode[STIX]{x1D64B}_{BE}-{\textstyle \frac{1}{5}}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D711}_{1}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(1)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}h_{2}\unicode[STIX]{x1D64B}_{BE},\end{eqnarray}$$

where the gradient operator $\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(\unicode[STIX]{x1D6FC})}$ acts only on velocity statistics of type- $\unicode[STIX]{x1D6FC}$ particles (e.g. $\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(1)}\boldsymbol{u}_{1}=\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{u}_{1}$ but $\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(1)}\boldsymbol{u}_{2}=\mathbf{0}$ ). Note that $\unicode[STIX]{x1D64B}_{BE}$ involves only velocity statistics (and not $\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}$ ). Thus, for example, we have the identity

(3.9) $$\begin{eqnarray}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(2)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}h_{1}\unicode[STIX]{x1D64B}_{BE}=\unicode[STIX]{x1D71A}_{1}h_{1}\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{BE}+\unicode[STIX]{x1D64B}_{BE}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{x}}\unicode[STIX]{x1D71A}_{1}h_{1},\end{eqnarray}$$

which can be used to evaluate the terms in $\boldsymbol{C}$ when required. Similar identities hold for the BE tensors in the energy balances below.

Comparing the momentum balances in (3.4) to those derived by Zhang et al. (Reference Zhang, Ma and Rauenzahn2006), we see that the interaction force densities, denoted by $n^{(i)}\overline{\boldsymbol{f}^{(ij)}}$ , have the following definitions for hard-sphere collisions:

(3.10) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}n^{(1)}\overline{\boldsymbol{f}^{(11)}}=\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D711}_{1}\unicode[STIX]{x1D748}^{(11)}=-\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{1},\quad n^{(1)}\overline{\boldsymbol{f}^{(12)}}=-\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{12}+\boldsymbol{C}-K\boldsymbol{u}_{12},\\ ~n^{(2)}\overline{\boldsymbol{f}^{(22)}}=\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D711}_{2}\unicode[STIX]{x1D748}^{(22)}=-\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{2},\quad n^{(2)}\overline{\boldsymbol{f}^{(21)}}=-\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{21}-\boldsymbol{C}+K\boldsymbol{u}_{12}.\end{array}\right\}\end{eqnarray}$$

The intraphase stress tensors $\unicode[STIX]{x1D748}^{(ii)}$ are thus related to the like-particle pressure tensors $\unicode[STIX]{x1D64B}_{i}$ , while the interphase stress tensors can be expressed as $\unicode[STIX]{x1D748}^{(ij)}=-\unicode[STIX]{x1D711}_{j}(\unicode[STIX]{x1D70C}_{2}\unicode[STIX]{x1D707}_{1}+\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D707}_{2})\frac{1}{5}{\mathcal{V}}_{12}g_{12}\unicode[STIX]{x1D64B}_{BE}$ where $\unicode[STIX]{x1D711}_{i}\unicode[STIX]{x1D748}^{(ij)}=\unicode[STIX]{x1D711}_{j}\unicode[STIX]{x1D748}^{(ji)}$ (Zhang et al. Reference Zhang, Ma and Rauenzahn2006). As discussed in Zhang et al. (Reference Zhang, Ma and Rauenzahn2006), the sum of the interspecies force densities is not equal to zero because they are calculated for different particle configurations, thereby giving rise to the interspecies stress tensor. For hard-sphere collisions of finite-size particles, this stress is proportional to $\unicode[STIX]{x1D64B}_{BE}$ .

The momentum-exchange term in (3.8) produces the Archimedes force of phase 1 on phase 2 when $d_{1}\ll d_{2}$ (i.e. $\unicode[STIX]{x1D707}_{1}=0$ ), as well as other contributions due to density and velocity gradients. As shown by Barajas et al. (Reference Barajas, Garcia-Colin and Pinã1973), the mixture pressure contribution from unlike collisions ( $\unicode[STIX]{x1D64B}_{12}+\unicode[STIX]{x1D64B}_{21}$ ) does not depend on $y$ . However, as seen above, the phase pressure tensors do depend on the choice of  $y$ . Finally, it is remarkable that $\unicode[STIX]{x1D64B}_{BE}$ is always anisotropic due to the mean-slip velocity. Momentum-exchange terms involving gradients of the mean-slip-velocity kinetic energy are known to result in a hyperbolic system under certain conditions (see Syamlal Reference Syamlal2011, in the context of a simpler two-fluid model).

3.3 Energy tensor balances

Neglecting the third-order central moments $\unicode[STIX]{x1D64C}$ , the balance equations for the total energy tensors $\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D640}_{\unicode[STIX]{x1D6FC}}$ found from (2.1) are

(3.11) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D640}_{1}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\boldsymbol{u}_{1}\otimes \unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D640}_{1}+(\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}\otimes \boldsymbol{u}_{1})^{\ast }+\unicode[STIX]{x1D642}_{1}+\unicode[STIX]{x1D642}_{12})=\boldsymbol{D}_{1}+\boldsymbol{D}_{0,12}+\boldsymbol{D}_{1,12},\\ \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D640}_{2}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\boldsymbol{u}_{2}\otimes \unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D640}_{2}+(\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D748}_{2}\otimes \boldsymbol{u}_{2})^{\ast }+\unicode[STIX]{x1D642}_{2}+\unicode[STIX]{x1D642}_{21})=\boldsymbol{D}_{2}+\boldsymbol{D}_{0,21}+\boldsymbol{D}_{1,21},\end{array}\right\}\end{eqnarray}$$

where $X_{ijk}^{\ast }:=\frac{1}{2}(\unicode[STIX]{x1D61F}_{ijk}+\unicode[STIX]{x1D61F}_{ikj})$ for the third-order tensor $\unicode[STIX]{x1D653}$ with component $\unicode[STIX]{x1D61F}_{ijk}$ . On the right-hand side of (3.11), the traceless isotropization terms, which result directly from the unclosed collision integral due to like collisions, are defined by

(3.12) $$\begin{eqnarray}\boldsymbol{D}_{\unicode[STIX]{x1D6FC}}:={\displaystyle \frac{6g_{0}\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FC}}\sqrt{\unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D6FC}}}}{\sqrt{\unicode[STIX]{x03C0}}d_{\unicode[STIX]{x1D6FC}}}}(\unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D748}_{\unicode[STIX]{x1D6FC}}).\end{eqnarray}$$

In the limit $d_{1}\ll d_{2}$ , $\boldsymbol{D}_{1}$ will drive $\unicode[STIX]{x1D748}_{1}$ towards isotropy despite the other anisotropic source terms, in which case it suffices to solve for the scalar energy $E_{1}$ .

The exchange terms in (3.11) due to unlike collisions have two parts. The first part does not involve spatial gradients

(3.13) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{D}_{0,12}:=K[\frac{1}{4}\unicode[STIX]{x1D707}_{2}(\unicode[STIX]{x1D6E9}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D748}_{2}+\boldsymbol{u}_{12}\otimes \boldsymbol{u}_{12})-\unicode[STIX]{x1D748}_{1}-(\boldsymbol{u}_{1}\otimes \boldsymbol{u}_{12})^{\dagger }],\\ \boldsymbol{D}_{0,21}:=K[\frac{1}{4}\unicode[STIX]{x1D707}_{1}(\unicode[STIX]{x1D6E9}\unicode[STIX]{x1D644}+\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D748}_{2}+\boldsymbol{u}_{12}\otimes \boldsymbol{u}_{12})-\unicode[STIX]{x1D748}_{2}+(\boldsymbol{u}_{2}\otimes \boldsymbol{u}_{12})^{\dagger }],\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D653}^{\dagger }:=\frac{1}{2}(\unicode[STIX]{x1D653}+\unicode[STIX]{x1D653}^{\text{t}})$ is the symmetric second-order tensor formed from $\unicode[STIX]{x1D653}$ . Note that $\boldsymbol{D}_{0,12}+\boldsymbol{D}_{0,21}=\frac{1}{2}K(\unicode[STIX]{x1D6E9}\unicode[STIX]{x1D644}-\unicode[STIX]{x1D748}_{1}-\unicode[STIX]{x1D748}_{2}-\boldsymbol{u}_{12}\otimes \boldsymbol{u}_{12})$ has zero trace, and thus these terms do not alter the total energy of the system. Furthermore, they do not directly affect the hyperbolic nature of the two-fluid model.

As shown in appendix A, the third-order collisional flux tensors in (3.11) have the forms

(3.14) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D642}_{1}:=(\unicode[STIX]{x1D64B}_{1}\otimes \boldsymbol{u}_{1})^{\ast },\quad \unicode[STIX]{x1D642}_{2}:=(\unicode[STIX]{x1D64B}_{2}\otimes \boldsymbol{u}_{2})^{\ast },\\ \unicode[STIX]{x1D642}_{12}:=\frac{1}{5}{\mathcal{V}}_{12}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D711}_{2}g_{12}\unicode[STIX]{x1D642}_{BE,12},\quad \unicode[STIX]{x1D642}_{21}:=\frac{1}{5}{\mathcal{V}}_{12}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D711}_{1}g_{12}\unicode[STIX]{x1D642}_{BE,21},\end{array}\right\}\end{eqnarray}$$

where the third-order tensors $\unicode[STIX]{x1D642}_{BE}$ due to the BE contribution (given in table 1) are expressed exactly in terms of the velocity moments. As expected, the components of $\unicode[STIX]{x1D642}_{BE}$ have the property $\unicode[STIX]{x1D60E}_{BE,ijk}=\unicode[STIX]{x1D60E}_{BE,ikj}$ . In terms of these tensors, the gradient-dependent contributions to the energy-exchange terms are

(3.15) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{D}_{1,12}:=\frac{1}{5}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D711}_{2}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(2)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}h_{1}\unicode[STIX]{x1D642}_{BE,12}-\frac{1}{5}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D711}_{1}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(1)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}h_{2}\unicode[STIX]{x1D642}_{BE,12},\\ ~\boldsymbol{D}_{1,21}:=\frac{1}{5}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D711}_{2}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(2)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}h_{1}\unicode[STIX]{x1D642}_{BE,21}-\frac{1}{5}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D711}_{1}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(1)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}h_{2}\unicode[STIX]{x1D642}_{BE,21}.\end{array}\right\}\end{eqnarray}$$

Using the definition of $\unicode[STIX]{x1D642}_{BE}$ , it can be shown that $\text{tr}(\boldsymbol{D}_{1,12})+\text{tr}(\boldsymbol{D}_{1,21})=0$ so that the total energy of the system is conserved.

Table 1.

Third-order tensors $\unicode[STIX]{x1D642}_{BE}$ and contracted first-order tensor $\unicode[STIX]{x1D666}_{BE}$ appearing in the energy conservation equations. In these expressions, the vectors $\boldsymbol{u}_{1}$ , $\boldsymbol{u}_{2}$ , $\boldsymbol{u}_{12}$ and $\boldsymbol{u}_{21}=-\boldsymbol{u}_{12}$ , and the velocity covariance tensors $\unicode[STIX]{x1D748}_{1}$ and $\unicode[STIX]{x1D748}_{2}$ appear. For isotropic systems, $\unicode[STIX]{x1D748}_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D644}$ and the formulae simplify considerably. The third-order tensors have the symmetry property $\unicode[STIX]{x1D60E}_{ijk}=\unicode[STIX]{x1D60E}_{ikj}$ . Using these definitions and $\unicode[STIX]{x1D707}_{1}+\unicode[STIX]{x1D707}_{2}=2$ , it is straightforward to prove that $\unicode[STIX]{x1D666}_{BE,12}=\unicode[STIX]{x1D666}_{BE,21}$ . Thus, we define $\unicode[STIX]{x1D666}_{BE}:=\frac{1}{2}(\unicode[STIX]{x1D666}_{BE,12}+\unicode[STIX]{x1D666}_{BE,21})$ . As a corollary, $\unicode[STIX]{x1D642}_{BE,12}-\unicode[STIX]{x1D642}_{BE,21}=\frac{1}{14}[8\langle \boldsymbol{g}\otimes \boldsymbol{g}\otimes \boldsymbol{g}\rangle +(\unicode[STIX]{x1D644}\otimes \langle g^{2}\boldsymbol{g}\rangle )^{\ast }-3\langle g^{2}\boldsymbol{g}\rangle \otimes \unicode[STIX]{x1D644}]$ and contraction of the difference is null. Finally, note that $\langle g^{2}\boldsymbol{g}\rangle +|\boldsymbol{u}_{12}|^{2}\boldsymbol{u}_{21}$ is proportional to $\unicode[STIX]{x1D64B}_{BE}\boldsymbol{\cdot }\boldsymbol{u}_{12}$ , which can be used to re-express $\unicode[STIX]{x1D666}_{BE}$ with terms involving $\unicode[STIX]{x1D64B}_{BE}$ .

3.4 Total energy balances

Balances for the total energy in each phase, $\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}E_{\unicode[STIX]{x1D6FC}}:=\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}\text{tr}(\unicode[STIX]{x1D640}_{\unicode[STIX]{x1D6FC}})$ , can be found by taking the trace of (3.11). This leads to

(3.16) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{1}E_{1}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}E_{1}+\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}\boldsymbol{\cdot }\boldsymbol{u}_{1}+\unicode[STIX]{x1D64B}_{1}\boldsymbol{\cdot }\boldsymbol{u}_{1}+\unicode[STIX]{x1D666}_{12})=C_{E}+D_{E},\\ ~\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{2}E_{2}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}E_{2}+\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D748}_{2}\boldsymbol{\cdot }\boldsymbol{u}_{2}+\unicode[STIX]{x1D64B}_{2}\boldsymbol{\cdot }\boldsymbol{u}_{2}+\unicode[STIX]{x1D666}_{21})=-C_{E}-D_{E},\end{array}\right\}\end{eqnarray}$$

where the drag exchange term is given by (as done for $\unicode[STIX]{x1D666}_{BE}$ in table 1, we use the symmetric form for $D_{E}$ )

(3.17) $$\begin{eqnarray}D_{E}:={\textstyle \frac{1}{2}}(\text{tr}(\boldsymbol{D}_{0,12})-\text{tr}(\boldsymbol{D}_{0,21}))={\textstyle \frac{3}{2}}K[\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D6E9}_{2}-\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D6E9}_{1}-{\textstyle \frac{1}{3}}(\unicode[STIX]{x1D707}_{2}\boldsymbol{u}_{2}+\unicode[STIX]{x1D707}_{1}\boldsymbol{u}_{1})\boldsymbol{\cdot }\boldsymbol{u}_{12}],\end{eqnarray}$$

and the energy fluxes by

(3.18a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D666}_{12}:={\textstyle \frac{1}{5}}{\mathcal{V}}_{12}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D711}_{2}g_{12}\unicode[STIX]{x1D666}_{BE},\quad \unicode[STIX]{x1D666}_{21}:={\textstyle \frac{1}{5}}{\mathcal{V}}_{12}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D711}_{1}g_{12}\unicode[STIX]{x1D666}_{BE}.\end{eqnarray}$$

The BE collisional correction to the energy flux $\unicode[STIX]{x1D666}_{BE}$ has components $\unicode[STIX]{x1D632}_{BE,i}:=\sum _{j}\unicode[STIX]{x1D60E}_{BE,ijj}$ (see table 1). The total-energy-exchange term is given by

(3.19) $$\begin{eqnarray}C_{E}:=\text{tr}(\boldsymbol{D}_{1,12})={\textstyle \frac{1}{5}}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D711}_{2}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(2)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}h_{1}\unicode[STIX]{x1D666}_{BE}-{\textstyle \frac{1}{5}}\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D711}_{1}(\unicode[STIX]{x2202}_{\boldsymbol{x}}-\unicode[STIX]{x2202}_{\boldsymbol{x}}^{(1)})\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}h_{2}\unicode[STIX]{x1D666}_{BE}.\end{eqnarray}$$

Thus, we can observe that the system total energy $\unicode[STIX]{x1D71A}E:=\unicode[STIX]{x1D71A}_{1}E_{1}+\unicode[STIX]{x1D71A}_{2}E_{2}$ is a conserved quantity. It is worth noting the important contribution of the relative velocity $\boldsymbol{u}_{12}$ in the production and exchange of energy between the two phases. At equilibrium, the mass ratios $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D6FC}}$ determine the energy distribution between the two phases.

3.5 Covariance tensors and fluctuating energy

Balance equations for the covariance tensors can be derived by subtracting off the contributions from the mean kinetic energy tensor $\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D646}_{\unicode[STIX]{x1D6FC}}:=\frac{1}{2}\unicode[STIX]{x1D71A}_{\unicode[STIX]{x1D6FC}}\boldsymbol{u}_{\unicode[STIX]{x1D6FC}}\otimes \boldsymbol{u}_{\unicode[STIX]{x1D6FC}}$ in (3.11). First, starting from the mean momentum balances, we derive balances for the mean kinetic energy tensors

(3.20) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D646}_{1}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}\otimes \unicode[STIX]{x1D646}_{1}+(\boldsymbol{u}_{1}\otimes \unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D64B}_{1}+\unicode[STIX]{x1D64B}_{12}))^{\dagger }~\\ \quad =(\boldsymbol{u}_{1}\otimes \boldsymbol{C}-K\boldsymbol{u}_{1}\otimes \boldsymbol{u}_{12})^{\dagger },~\\ \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{p}\unicode[STIX]{x1D646}_{2}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}\otimes \unicode[STIX]{x1D646}_{2}+(\boldsymbol{u}_{2}\otimes \unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }(\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D748}_{2}+\unicode[STIX]{x1D64B}_{2}+\unicode[STIX]{x1D64B}_{21}))^{\dagger }\\ ~\quad =-(\boldsymbol{u}_{2}\otimes \boldsymbol{C}-K\boldsymbol{u}_{2}\otimes \boldsymbol{u}_{12})^{\dagger }.\end{array}\right\}\end{eqnarray}$$

We then subtract these expressions from the total energy tensor balances to find

(3.21) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}+{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}\otimes \unicode[STIX]{x1D748}_{1}+((\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D64B}_{1})\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{u}_{1})^{\dagger }+\unicode[STIX]{x2202}_{\boldsymbol{ x}}\boldsymbol{\cdot }\unicode[STIX]{x1D642}_{12}-(\boldsymbol{u}_{1}\otimes \unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{12})^{\dagger }\\ ~\quad =\boldsymbol{D}_{1}+\boldsymbol{D}_{0,12}+K(\boldsymbol{u}_{1}\otimes \boldsymbol{u}_{12})^{\dagger }+\boldsymbol{D}_{1,12}-(\boldsymbol{u}_{1}\otimes \boldsymbol{C})^{\dagger },\\ ~{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D748}_{2}+{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}\otimes \unicode[STIX]{x1D748}_{2}+((\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D748}_{2}+\unicode[STIX]{x1D64B}_{2})\boldsymbol{\cdot }\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{u}_{2})^{\dagger }+\unicode[STIX]{x2202}_{\boldsymbol{ x}}\boldsymbol{\cdot }\unicode[STIX]{x1D642}_{21}-(\boldsymbol{u}_{2}\otimes \unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{21})^{\dagger }\\ ~\quad =\boldsymbol{D}_{2}+\boldsymbol{D}_{0,21}-K(\boldsymbol{u}_{2}\otimes \boldsymbol{u}_{12})^{\dagger }+\boldsymbol{D}_{1,21}+(\boldsymbol{u}_{2}\otimes \boldsymbol{C})^{\dagger }.\end{array}\right\}\end{eqnarray}$$

The velocity covariance tensors will be non-zero if there is a non-zero mean-slip velocity. Also, in general, they will be anisotropic due to the anisotropy of their source terms.

The traces of the velocity covariance tensors in (3.21) lead to the fluctuating energy balance equations

(3.22) $$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\frac{3}{2}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D6E9}_{1}+\frac{3}{2}\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}\unicode[STIX]{x1D6E9}_{1}+(\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D64B}_{1})\boldsymbol{ : }\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{u}_{1}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D666}_{12}-\boldsymbol{u}_{1}\boldsymbol{\cdot }(\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{12})\\[3.0pt] \quad =D_{E}+K\boldsymbol{u}_{1}\boldsymbol{\cdot }\boldsymbol{u}_{12}+C_{E}-\boldsymbol{u}_{1}\boldsymbol{\cdot }\boldsymbol{C},\\[3.0pt] \frac{3}{2}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D6E9}_{2}+\frac{3}{2}\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}\unicode[STIX]{x1D6E9}_{2}+(\unicode[STIX]{x1D71A}_{2}\unicode[STIX]{x1D748}_{2}+\unicode[STIX]{x1D64B}_{2})\boldsymbol{ : }\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{u}_{2}+\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D666}_{21}-\boldsymbol{u}_{2}\boldsymbol{\cdot }(\unicode[STIX]{x2202}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{21})\\[3.0pt] \quad =-D_{E}-K\boldsymbol{u}_{2}\boldsymbol{\cdot }\boldsymbol{u}_{12}-C_{E}+\boldsymbol{u}_{2}\boldsymbol{\cdot }\boldsymbol{C},\end{array}\right\}\end{eqnarray}$$

which govern the phase temperatures. The terms on the left-hand sides involving the phase pressure tensors correspond to compression heating. The other terms on the right-hand side are production/redistribution of fluctuating energy between phases. The sum of the two balances in (3.22) is the total fluctuating energy. For the latter, we can observe that there is a source term $K|\boldsymbol{u}_{12}|^{2}$ due to the transfer of mean kinetic energy to fluctuations when the phases have a mean-slip velocity. On the other hand, $\boldsymbol{u}_{21}\boldsymbol{\cdot }\boldsymbol{C}$ can function as a source or sink of fluctuating energy depending on the signs of the gradients appearing in $\boldsymbol{C}$ . These observations imply that even in the absence of mean gradients, the phase temperatures will be non-zero if the phases are initially in relative mean motion. Finally, using the momentum balances in (3.4) and (3.22), it is possible to derive the balance equation for the ‘mixture’ temperature $\unicode[STIX]{x1D6E9}$ .

3.6 Equilibrium state

The equilibrium state for a homogeneous, binary system is $\boldsymbol{u}_{1}=\boldsymbol{u}_{2}$ , $\unicode[STIX]{x1D707}_{1}\unicode[STIX]{x1D6E9}_{1}=\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D6E9}_{2}$ and $\unicode[STIX]{x1D748}_{\unicode[STIX]{x1D6FC}}=\unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D6FC}}\unicode[STIX]{x1D644}$ (see, e.g. Kolodner Reference Kolodner1957; Grad Reference Grad and Devienne1960; Morse Reference Morse1963; Goldman & Sirovich Reference Goldman and Sirovich1967), and the densities, velocities and temperatures are spatially uniform. The functions $h_{1}$ and $h_{2}$ introduced in (2.6) are chosen such that the isotropic flux pressure, i.e. for type-1 particles

(3.23) $$\begin{eqnarray}\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D748}_{1}+\unicode[STIX]{x1D64B}_{1}+\unicode[STIX]{x1D64B}_{12}=\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D6E9}_{1}(1+{\mathcal{V}}_{11}\unicode[STIX]{x1D711}_{1}g_{0}+{\mathcal{V}}_{12}\unicode[STIX]{x1D711}_{2}g_{12})\unicode[STIX]{x1D644}=p_{1}\unicode[STIX]{x1D644},\end{eqnarray}$$

and the buoyancy pressure, defined to be the non-zero term in $\boldsymbol{C}$ at equilibrium, i.e. for type-1 particles (the other terms in $\boldsymbol{C}$ depend on gradients in density and velocity and correspond to forces other than the buoyancy force)

(3.24) $$\begin{eqnarray}{\textstyle \frac{1}{5}}\unicode[STIX]{x1D707}_{2}\unicode[STIX]{x1D71A}_{1}h_{1}\unicode[STIX]{x1D64B}_{BE}|_{eq}=\unicode[STIX]{x1D71A}_{1}\unicode[STIX]{x1D6E9}_{1}h_{1}\unicode[STIX]{x1D644}=p_{1}\unicode[STIX]{x1D644},\end{eqnarray}$$

are identical for each particle type. In equilibrium, the phase pressures are isotropic and the momentum balances reduce to

(3.25a,b ) $$\begin{eqnarray}\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{1}=\unicode[STIX]{x1D711}_{2}\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{1}-\unicode[STIX]{x1D711}_{1}\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{2},\quad \unicode[STIX]{x2202}_{\boldsymbol{x}}p_{2}=\unicode[STIX]{x1D711}_{1}\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{2}-\unicode[STIX]{x1D711}_{2}\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{1}.\end{eqnarray}$$

Thus, $\unicode[STIX]{x2202}_{\boldsymbol{x}}(p_{1}+p_{2})=\mathbf{0}$ , i.e. the mixture pressure is uniform, and hence $\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{1}=-\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{2}$ . It then follows that in order for the above momentum balances to hold, $\unicode[STIX]{x2202}_{\boldsymbol{x}}p_{\unicode[STIX]{x1D6FC}}=\mathbf{0}$ at equilibrium, i.e. the phase pressures must be uniform.

If the mean momentum is initially null, $\unicode[STIX]{x1D71A}_{1}\boldsymbol{u}_{1}+\unicode[STIX]{x1D71A}_{2}\boldsymbol{u}_{2}=\mathbf{0}$ but $\boldsymbol{u}_{1}\neq \boldsymbol{u}_{2}$ , then at equilibrium $\boldsymbol{u}_{1}=\boldsymbol{u}_{2}=\mathbf{0}$ and the mean kinetic energy will be transferred to $\unicode[STIX]{x1D6E9}_{\unicode[STIX]{x1D6FC}}$ through collisions. Thus, in the limit where $\unicode[STIX]{x1D707}_{1}\ll \unicode[STIX]{x1D707}_{2}$ , the phase temperatures will be very different at equilibrium with $\unicode[STIX]{x1D6E9}_{2}\ll \unicode[STIX]{x1D6E9}_{1}$ . (For systems in motion, $\unicode[STIX]{x1D6E9}_{2}$ has a source term due to phase compression and, hence, can be much larger than its equilibrium value. For fluid–particle systems in § 4, the equilibrium value of $\unicode[STIX]{x1D6E9}_{2}$ is essentially zero.) As the particles are assumed to have no internal degrees of freedom, their internal energies are null. In the absence of external forces (e.g. gravity), the mean momentum and total energy of the system are conserved. As first described by Grad (Reference Grad and Devienne1960), the rates of relaxation to equilibrium for the mean velocities and granular temperatures can be very different for components with disparate masses. For example, under certain conditions (e.g. plasmas) the mean velocities will rapidly equilibrate so that only the two temperatures differ (i.e. the so-called two-temperature model (Sirovich Reference Sirovich1962)). Hereinafter, we will focus on cases for which a two-fluid description is required, i.e. the heavy component is a trace component with very small molar concentration as first described in Hamel (Reference Hamel1966).

5.1 Uncoupled model

Formally, the 1-D model can be rewritten as

(5.1) $$\begin{eqnarray}\unicode[STIX]{x1D63C}(\unicode[STIX]{x1D653})\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D653}+\unicode[STIX]{x1D63D}(\unicode[STIX]{x1D653})\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D653}=\boldsymbol{C}(\unicode[STIX]{x1D653}),\end{eqnarray}$$

where the coefficient matrices $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}:=\unicode[STIX]{x1D63D}^{\ast }+\unicode[STIX]{x1D63D}_{0}$ yield the flux matrix $\unicode[STIX]{x1D641}:=\unicode[STIX]{x1D63C}^{-1}\unicode[STIX]{x1D63D}$ . The matrix $\unicode[STIX]{x1D63D}_{0}$ is the contribution due to the coupled behaviour, and its form is complicated due to the nonlinear dependencies of $g_{0}$ and $g_{12}$ on $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FC}}$ . However, as done here, $\unicode[STIX]{x1D63D}_{0}$ can be easily computed using symbolic software. $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}^{\ast }$ are from the uncoupled system with only drag (i.e. $K>0$ and $P^{\ast }=q^{\ast }=g_{0}=0$ ). Written in terms of the components of $\unicode[STIX]{x1D653}$ , the latter are

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D63C}:=\left[\begin{array}{@{}cccccc@{}}1 & 0 & 0 & 0 & 0 & 0\\ 0 & \unicode[STIX]{x1D711}_{1} & 0 & 0 & 0 & 0\\ 0 & 0 & \unicode[STIX]{x1D711}_{1} & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & \unicode[STIX]{x1D711}_{2} & 0\\ 0 & 0 & 0 & 0 & 0 & \unicode[STIX]{x1D711}_{2}\end{array}\right],\end{eqnarray}$$

and

(5.3) $$\begin{eqnarray}\unicode[STIX]{x1D63D}^{\ast }:=\left[\begin{array}{@{}cc@{}}\unicode[STIX]{x1D63D}_{1} & 0\\ 0 & \unicode[STIX]{x1D63D}_{2}\end{array}\right],\end{eqnarray}$$

with

(5.4) $$\begin{eqnarray}\unicode[STIX]{x1D63D}_{1}:=\left[\begin{array}{@{}ccc@{}}u_{1} & \unicode[STIX]{x1D711}_{1} & 0\\[4.0pt] \frac{2}{3}E_{1}-\frac{1}{3}u_{1}^{2} & \frac{1}{3}\unicode[STIX]{x1D711}_{1}u_{1} & \frac{2}{3}\unicode[STIX]{x1D711}_{1}\\[4.0pt] u_{1}\left(\frac{2}{3}E_{1}-\frac{1}{3}u_{1}^{2}\right) & \unicode[STIX]{x1D711}_{1}\left(\frac{2}{3}E_{1}-u_{1}^{2}\right) & \frac{5}{3}\unicode[STIX]{x1D711}_{1}u_{1}\end{array}\right],\end{eqnarray}$$

and

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D63D}_{2}:=\left[\begin{array}{@{}ccc@{}}u_{2} & \unicode[STIX]{x1D711}_{2} & 0\\[4.0pt] \frac{2}{3}E_{2}-\frac{1}{3}u_{2}^{2} & \frac{1}{3}\unicode[STIX]{x1D711}_{2}u_{2} & \frac{2}{3}\unicode[STIX]{x1D711}_{2}\\[4.0pt] u_{2}\left(\frac{2}{3}E_{2}-\frac{1}{3}u_{2}^{2}\right) & \unicode[STIX]{x1D711}_{2}\left(\frac{2}{3}E_{2}-u_{2}^{2}\right) & \frac{5}{3}\unicode[STIX]{x1D711}_{2}u_{2}\end{array}\right].\end{eqnarray}$$

Given the form of (5.3), we observe that the fluid- and particle-phase fluxes are uncoupled and have an identical structure. The six eigenvalues for the flux matrix of the uncoupled model are

(5.6a,b ) $$\begin{eqnarray}u_{1},u_{1}\pm \sqrt{{\textstyle \frac{5}{3}}\unicode[STIX]{x1D6E9}_{1}},\quad u_{2},u_{2}\pm \sqrt{{\textstyle \frac{5}{3}}\unicode[STIX]{x1D6E9}_{2}},\end{eqnarray}$$

and hence the uncoupled model is hyperbolic if $0<\unicode[STIX]{x1D6E9}_{1}$ and $0<\unicode[STIX]{x1D6E9}_{2}$ . Given that the mass ratio between the particles and fluid molecules is very large, we expect to have $\unicode[STIX]{x1D6E9}_{2}\ll \unicode[STIX]{x1D6E9}_{1}$ (i.e. the speed of sound in the fluid phase is much larger than in the particle phase). As expected (Marble Reference Marble1970), with only drag, the flux matrix for a very dilute system (i.e. $\unicode[STIX]{x1D711}\ll 1$ ) has the form of two uncoupled compressible (i.e. perfect gas) systems with independent eigenvalues. Recall, however, that when $\unicode[STIX]{x1D711}_{2}\ll 1$ , particle–particle collisions will be rare, and therefore second-order models for the particle phase will be inaccurate unless the particle Stokes number is very small (see discussion in Vié et al. Reference Vié, Doisneau and Massot2015).

5.2 Application of Sturm’s theorem

The eigenvalues of the complete 1-D model depend on the following five parameters:

(5.7a-e ) $$\begin{eqnarray}Ma_{s}={\displaystyle \frac{u_{2}-u_{1}}{\sqrt{5\unicode[STIX]{x1D6E9}_{1}/3}}},\quad \unicode[STIX]{x1D6E9}_{r}={\displaystyle \frac{\unicode[STIX]{x1D6E9}_{2}}{\unicode[STIX]{x1D6E9}_{1}}},\quad {\mathcal{Z}}={\displaystyle \frac{\unicode[STIX]{x1D70C}_{f}}{\unicode[STIX]{x1D70C}_{p}}}={\displaystyle \frac{\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D711}_{1}}{\unicode[STIX]{x1D70C}_{2}(1-\unicode[STIX]{x1D711}_{2})}},\quad \unicode[STIX]{x1D711}_{1},\quad \unicode[STIX]{x1D711}_{2}.\end{eqnarray}$$

Because the model is Galilean invariant, in the frame of reference of the fluid, we can choose $u_{1}=0$ in the characteristic polynomial for the eigenvalues. Using this frame of reference, the six eigenvalues of the uncoupled model are

(5.8) $$\begin{eqnarray}0,\pm 1,Ma_{s},Ma_{s}\pm \sqrt{\unicode[STIX]{x1D6E9}_{r}},\end{eqnarray}$$

and will be used to compare with the eigenvalues of the complete model. For the latter, the flux matrix $\unicode[STIX]{x1D641}$ and the characteristic polynomial $p_{0}(x)$ can be computed using symbolic software, and depend on the five parameters above. In order to determine whether this sixth-degree polynomial has real roots, we make use of Sturm’s theorem (Sturm Reference Sturm1829) for the interval $(-\infty ,+\infty )$ . For this purpose, we compute symbolically the Sturm chain $p_{0}(x),p_{1}(x),\ldots ,p_{6}(x)$ using polynomial long division (https://fr.wikipedia.org/wiki/Théorème_de_Sturm). This sequence of polynomials has decreasing degree such that $p_{6}$ is zero-degree. According to Sturm’s theorem, it suffices that the coefficient of the highest-order term in each of these polynomials be positive to prove that all six eigenvalues are real valued. For the 1-D model, it is found that only $p_{5}$ and $p_{6}$ are needed to determine whether the system is hyperbolic for a given choice of ( $Ma_{s},\unicode[STIX]{x1D6E9}_{r},{\mathcal{Z}},\unicode[STIX]{x1D711}_{1},\unicode[STIX]{x1D711}_{2}$ ). While these two coefficients can be computed analytically using symbolic software, they are too complicated to analyse directly to find the constraints on the five parameters. However, they are easily computed numerically for a given set of parameters. In the analysis below, we define a test function $f(\unicode[STIX]{x1D711}_{2},\unicode[STIX]{x1D6E9}_{r})$ equal to the minimum of the Sturm coefficients for fixed parameters ( $Ma_{s},{\mathcal{Z}},\unicode[STIX]{x1D711}_{1}$ ). According to Sturm’s theorem, if $f(\unicode[STIX]{x1D711}_{2},\unicode[STIX]{x1D6E9}_{r})<0$ , then at least two eigenvalues are complex. Note that the maximum of $f$ is unity due to the normalization of $p_{0}$ . Below, we will use hyperbolicity plots in the $\unicode[STIX]{x1D711}_{2}$ – $\unicode[STIX]{x1D6E9}_{r}$ plane to show where the Sturm test function is non-negative for fixed values of ( $\unicode[STIX]{x1D711}_{1},Ma_{s},{\mathcal{Z}}$ ).

In general, as will be shown in detail below, the hyperbolic region of the two-fluid model in table 4 is similar to that found by Houim & Oran (Reference Houim and Oran2016). In analogy to their notation, $Ma_{s}$ corresponds to $M_{slip}$ , $\unicode[STIX]{x1D6E9}_{r}$ to $C_{r}$ , and (roughly) ${\mathcal{Z}}$ to $z$ . However, the present model depends separately on ${\mathcal{Z}},\unicode[STIX]{x1D711}_{1},\unicode[STIX]{x1D711}_{2}$ , while, written in our notation, $z={\mathcal{Z}}\unicode[STIX]{x1D711}_{2}/(1-\unicode[STIX]{x1D711}_{2})$ and $\unicode[STIX]{x1D711}_{1}\rightarrow 0$ in Houim & Oran (Reference Houim and Oran2016), i.e. the two parameters are equal for a perfect gas when $\unicode[STIX]{x1D6FC}_{p}=\frac{1}{2}$ ; otherwise, in most physically relevant cases, $z<{\mathcal{Z}}$ . For $Ma_{s}=0$ or ${\mathcal{Z}}=0$ , both models are strictly hyperbolic. However, as shown below, in contrast to the model in Houim & Oran (Reference Houim and Oran2016), the two-fluid model in table 4 is hyperbolic for all values of the other parameters if $Ma_{s}$ is small enough, including ${\mathcal{Z}}\geqslant 0.1$ . For larger $Ma_{s}$ , the model is hyperbolic except for certain combinations of $\unicode[STIX]{x1D6E9}_{r}$ and $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D6FC}}$ . Note that in a gas–particle flow, we expect physically that ${\mathcal{Z}}\ll 1$ , i.e. the density of the gas will be (much) smaller than that of the solid. Having values of ${\mathcal{Z}}\geqslant 1$ can only occur when the fluid is a liquid (or with $\unicode[STIX]{x1D70C}_{1}\gg \unicode[STIX]{x1D70C}_{2}$ ), in which case $\unicode[STIX]{x1D711}_{1}\approx 0.63$ for a pure fluid. For a pure liquid, $\unicode[STIX]{x1D711}_{1}$ could be chosen to produce the correct speed of sound. For the hard-sphere fluid, the speed of sound increases due to $g_{0}$ as $\unicode[STIX]{x1D711}$ approaches unity. Nevertheless, a stiffened gas model, coupled with an added-mass model, may be a better choice for representing a compressible liquid. Thus, any complex eigenvalues observed when ${\mathcal{Z}}\geqslant 0.1$ and $\unicode[STIX]{x1D711}_{1}\ll 0.63$ are unphysical because of the unphysical choice of parameters. Moreover, when ${\mathcal{Z}}\geqslant 0.1$ , the drag force will quickly damp out $\unicode[STIX]{x1D6E9}_{2}$ so that only small values of $\unicode[STIX]{x1D6E9}_{r}$ are physically relevant. For practical applications, we will mainly be interested in the perfect gas model, which corresponds to the two-fluid model in table 3. Nonetheless, for completeness, we will consider all mathematically relevant values when testing the system for real eigenvalues.

Figure 1.

Hyperbolicity plot for perfect gas ( $\unicode[STIX]{x1D711}_{1}=0$ ) with ${\mathcal{Z}}=0.001$ and varying $Ma_{s}$ . The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues.

5.3 Perfect gas

For a perfect gas ( $\unicode[STIX]{x1D711}_{1}\rightarrow 0$ ), the regions where the Sturm test function is negative are shown in figures 1–4 for ${\mathcal{Z}}=0.001$ , 0.01 and 0.1, respectively. As noted above, the Sturm test function is non-negative for ${\mathcal{Z}}=0$ . In figure 1, we observe that the Sturm test function is non-negative for $M_{s}\leqslant 1/3$ for all values of $\unicode[STIX]{x1D711}_{2}$ . For $Ma_{s}\geqslant 2/3$ , a curve exists for which the Sturm test function is negative. In figure 1(c), this curve extends to lower values of $\unicode[STIX]{x1D719}_{2}$ . Similarly, in figure 1(d), a second curve is evident along which the 1-D model loses hyperbolicity. In order to understand the behaviour of the eigenvalues when the Sturm test function is negative, we plot them against $\unicode[STIX]{x1D6E9}_{r}$ in figure 2. In this figure, the ‘gas-phase’ eigenvalues correspond to ( $-1$ , 0, 1), while two of the ‘particle-phase’ eigenvalues grow in magnitude with increasing $\unicode[STIX]{x1D6E9}_{r}$ . The complex eigenvalues, when they exist, occur when one of these two eigenvalues overlap with the gas-phase eigenvalues. From figure 2, the complex eigenvalues occur at a single point (which is consistent with figure 1). Thus, for ${\mathcal{Z}}$ values typical of gas–particle flows we can expect that the loss of hyperbolicity will not be issue in numerical simulations, even if $M_{s}$ is larger than unity. Moreover, in applications to real particles, we expect $0<\unicode[STIX]{x1D6E9}_{r}\ll 1$ (see Tavenashad et al. Reference Tavenashad, Passalacqua, Fox and Subramaniam2019, for examples) and, hence, the particle-phase eigenvalue will only interact with the gas-phase eigenvalue located at 0.

Figure 2.

Real part of eigenvalues versus $\unicode[STIX]{x1D6E9}_{r}$ for perfect gas ( $\unicode[STIX]{x1D711}_{1}=0$ ) with ${\mathcal{Z}}=0.001$ . Values of $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ are chosen to cross curves with a negative Sturm test function in figure 1.

The results for larger ${\mathcal{Z}}$ are given in figures 3–5, where it can be seen that the non-hyperbolic regions grow larger as ${\mathcal{Z}}$ increases. For $Ma_{s}\leqslant 1/3$ , the system is hyperbolic for all physically relevant values of $\unicode[STIX]{x1D711}_{2}$ and $\unicode[STIX]{x1D6E9}_{f}$ . When $Ma_{s}\geqslant 2/3$ , the system can have a pair of complex eigenvalues in larger regions of $\unicode[STIX]{x1D711}_{2}$ – $\unicode[STIX]{x1D6E9}_{r}$ space, which are significantly larger for ${\mathcal{Z}}=0.1$ . Houim & Oran (Reference Houim and Oran2016) found similar behaviour for their model, but in their case complex eigenvalues were found for $\unicode[STIX]{x1D6E9}_{r}=0$ for all values of $Ma_{s}$ . Having complex eigenvalues with $\unicode[STIX]{x1D6E9}_{r}=0$ is particularly troublesome as this corresponds to the equilibrium value of $\unicode[STIX]{x1D6E9}_{2}$ . However, tests with a kinetic-based added-mass term (not shown here) indicate that the system recovers hyperbolicity for small $\unicode[STIX]{x1D6E9}_{r}$ because with added mass the three particle-phase eigenvalues are no longer equal at $\unicode[STIX]{x1D6E9}_{r}=0$ . In figure 5 the eigenvalues are plotted for parameter sets that exhibit non-hyperbolic behaviour. We can observe from figure 5 that the Sturm test function correctly predicts the values of $\unicode[STIX]{x1D6E9}_{r}$ where two eigenvalues have the same real parts (i.e. complex conjugates).

Figure 3.

Hyperbolicity plot for perfect gas ( $\unicode[STIX]{x1D711}_{1}=0$ ) with ${\mathcal{Z}}=0.01$ and varying $Ma_{s}$ . The Sturm test function is negative in black regions, indicating that 1-D model has complex eigenvalues.

Figure 4.

Hyperbolicity plot for perfect gas ( $\unicode[STIX]{x1D711}_{1}=0$ ) with ${\mathcal{Z}}=0.1$ and varying $Ma_{s}$ . The Sturm test function is negative in black regions, indicating that 1-D model has complex eigenvalues.

Figure 5.

Real part of eigenvalues versus $\unicode[STIX]{x1D6E9}_{r}$ for perfect gas ( $\unicode[STIX]{x1D711}_{1}=0$ ) with ${\mathcal{Z}}=0.01$ and 0.1. Values of $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ are chosen to cross curves with a negative Sturm test function in figures 3 and 4.

In summary, for a perfect gas the two-fluid model is strictly hyperbolic for all physically realistic values of the parameters if $Ma_{s}$ and ${\mathcal{Z}}$ are not too large. The value of $Ma_{s}\leqslant 1/3$ is not a rigorous upper bound for hyperbolicity; however, as discussed in the derivation of the model, one should not expect the pair distribution function $g_{12}$ in (2.11) to be accurate for larger values of $Ma_{s}$ . Thus, the non-hyperbolic regions seen in the figures above for $Ma_{s}>1/3$ should not be encountered in gas–particle simulations for which the proposed two-fluid model is physically applicable. In appendix C a simpler system with $g_{0}=0$ is investigated and confirms these conclusions. Indeed, with $g_{0}=0$ the hyperbolic region is much larger than when the particle phase has a non-ideal pressure model, which would seem counter-intuitive vis-a-vis the coupling with the gas phase.

5.4 Non-ideal gas

By setting $\unicode[STIX]{x1D711}=0.01$ , the gas becomes non-ideal and the gas-phase eigenvalues (with ${\mathcal{Z}}=0$ ) depend on $g_{0}$ and thus are slightly larger than for the perfect gas. The hyperbolicity plots found from the Sturm test functions are given in figures 6–8, and example eigenvalues in figure 9. In these figures, all other parameters are treated in the same manner as done for the perfect gas. Qualitatively, the non-ideal gas behaves very similarly to the perfect gas. The regions in parameter space where the system is hyperbolic are slightly modified by increasing $\unicode[STIX]{x1D711}_{1}$ . Quantitatively, the eigenvalues have larger magnitudes due to the higher gas-phase speed of sound. In addition, when $M_{s}>1/3$ there is a region with large $\unicode[STIX]{x1D711}_{2}$ and $\unicode[STIX]{x1D6E9}_{r}$ that has complex eigenvalues, and this region moves to smaller $\unicode[STIX]{x1D711}_{2}$ as ${\mathcal{Z}}$ increase. Overall, the two-fluid model with a non-ideal gas is again hyperbolic for physically realistic values of the parameters if $Ma_{s}$ and ${\mathcal{Z}}$ are not too large.

Figure 6.

Hyperbolicity plot for non-ideal gas ( $\unicode[STIX]{x1D711}_{1}=0.01$ ) with ${\mathcal{Z}}=0.001$ and varying $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ . The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues.

Figure 7.

Hyperbolicity plot for non-ideal gas ( $\unicode[STIX]{x1D711}_{1}=0.01$ ) with ${\mathcal{Z}}=0.01$ and varying $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ . The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues.

Figure 8.

Hyperbolicity plot for non-ideal gas ( $\unicode[STIX]{x1D711}_{1}=0.01$ ) with ${\mathcal{Z}}=0.1$ and varying $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ . The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues.

Figure 9.

Real part of eigenvalues versus $\unicode[STIX]{x1D6E9}_{r}$ for non-ideal gas ( $\unicode[STIX]{x1D711}_{1}=0.01$ ) with ${\mathcal{Z}}=0.001$ , 0.01, 0.1 and $Ma_{s}=1$ and $4/3$ . Values of $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ are chosen to cross curves with a negative Sturm test function in figures 6–8.

5.5 Dense gas

With $\unicode[STIX]{x1D711}_{1}=0.1$ , the gas is dense, i.e. the fluid molecular volume fraction is as large as the particle volume fraction $\unicode[STIX]{x1D711}_{2}$ . The Sturm test function plots for this case are given in figures 10–12, and example eigenvalues in figure 13. In these figures, all other parameters are treated in the same manner as done for the perfect and non-ideal gases. As expected, the effect of $g_{0}$ is to push the two ‘gas-phase’ eigenvalues further from the origin. From the hyperbolicity plots, we can observe that the dense gas generally makes the system less hyperbolic for a given value of $\unicode[STIX]{x1D711}_{2}$ . For example, comparing figure 10 to figure 1, we can observe that for $M_{s}>2/3$ the non-hyperbolic regions found with $\unicode[STIX]{x1D711}_{2}=0.1$ occur at larger $\unicode[STIX]{x1D6E9}_{r}$ in figure 10. Furthermore, in figure 12 we see that for the dense gas the non-hyperbolic region at $\unicode[STIX]{x1D6E9}_{r}=0$ is no longer present. In summary, the two-fluid model with a dense gas is hyperbolic for physically realistic values of the parameters if $Ma_{s}$ is not too large.

Figure 10.

Hyperbolicity plot for dense gas ( $\unicode[STIX]{x1D711}_{1}=0.1$ ) with ${\mathcal{Z}}=0.001$ and varying $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ . The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues.

Figure 11.

Hyperbolicity plot for dense gas ( $\unicode[STIX]{x1D711}_{1}=0.1$ ) with ${\mathcal{Z}}=0.01$ and varying $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ . The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues.

Figure 12.

Hyperbolicity plot for dense gas ( $\unicode[STIX]{x1D711}_{1}=0.1$ ) with ${\mathcal{Z}}=0.1$ and varying $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ . The Sturm test function is negative in black regions, indicating that the 1-D model has complex eigenvalues.

Figure 13.

Real part of eigenvalues versus $\unicode[STIX]{x1D6E9}_{r}$ for dense gas ( $\unicode[STIX]{x1D711}_{1}=0.1$ ) with ${\mathcal{Z}}=0.010$ , 0.01, 0.1 and $Ma_{s}=2/3$ , 1 and $4/3$ . Values of $Ma_{s}$ and $\unicode[STIX]{x1D711}_{2}$ are chosen to cross curves with a negative Sturm test function in figures 10–12.