1. Introduction
Despite the fact that such models are most often ill posed (Lhuillier, Chang & Theofanous Reference Lhuillier, Chang and Theofanous2013), two-fluid models have found widespread use for simulating industrial and environmental polydisperse multiphase flows (Fox Reference Fox2024). The principal modelling difficulty to achieving a well-posed system results from the Archimedes force, involving the continuous fluid pressure gradient, for cases where the particle-phase material density
$\rho _p$
is smaller than that of the fluid
$\rho _f$
(Lhuillier et al. Reference Lhuillier, Chang and Theofanous2013). In addition to being non-physical, ill-posed two-fluid models suffer from numerical instabilities due to complex eigenvalues, yielding spurious oscillations (Panicker, Passalacqua & Fox Reference Panicker, Passalacqua and Fox2018). Given that two-fluid models are usually derived from well-posed microscale models (e.g. the Navier–Stokes equation for the fluid phase coupled to Newton’s equations for the particles) (Drew & Passman Reference Drew and Passman1998; Zhang, Ma & Rauenzahn Reference Zhang, Ma and Rauenzahn2006), it is clear that they are ill posed due to the approximations used to close the model equations. In general, there are two types of unclosed terms: (i) terms modifying the spatial fluxes of mass, momentum and energy; and (ii) source terms for exchange of mass, momentum and energy. For example, in the volume-averaging approach (Drew & Passman Reference Drew and Passman1998), the source terms for momentum and energy exchange involve unclosed particle-surface integrals depending on the fluid stresses at the particle surface. Similarly, in the kinetic approach of Zhang et al. (Reference Zhang, Ma and Rauenzahn2006) an unclosed particle–fluid–particle (pfp) stress tensor
$\boldsymbol {\Sigma }_{pfp}$
arises from fluid-mediated (or hydrodynamic) forces on a given particle due to the surrounding particles (see Wang et al. Reference Wang, Yang, Zhang and Balachandar2021; Zhang Reference Zhang2021 for more details).
With these difficulties in mind, Fox (Reference Fox2019) derived a two-fluid model for ideal-fluid–particle flow starting from the Boltzmann–Enskog (BE) kinetic equation for hard spheres with disparate sizes. As done in the derivation of the Euler equations for gas dynamics (Cercignani Reference Cercignani1988), the goal was to find exact expressions for the fluxes and source terms in the two-fluid model for a specific system with well-defined interactions between the phases (i.e. binary, hard-sphere collisions). Such a model can then be analysed to check if it is well posed and, if so, to determine what are the necessary properties of the fluxes and source terms. As is well known, the (inviscid) Euler equations for an ideal fluid can be found from the Boltzmann kinetic equation for monodisperse hard spheres. For a binary system of hard spheres with different diameters
$d_1 \ll d_2$
, the BE kinetic equation is used to model how spatial gradients in the macroscopic variables (e.g. density, velocity, pressure) affect the collision rates between particles with different diameters (Chapman & Cowling Reference Chapman and Cowling1952). For example, in order to include the Archimedes force in the momentum balances, the binary collision operators in the BE kinetic equation require a specific form for the particle-pair distribution function (see Fox Reference Fox2019 and § 2 for details). In addition to the Archimedes force, this distribution function generates a momentum-exchange term (denoted by
$\boldsymbol {F}_{\!\!fp} = - \boldsymbol {F}_{\!pf}$
) involving spatial gradients in the fluid-phase density
$\rho _f$
and the fluid-phase velocity
$\boldsymbol{u}_f$
. Likewise, the presence of particles modifies the momentum and energy fluxes in the fluid phase. For the binary hard-sphere system, it is in fact these additional terms that are responsible for making the two-fluid model well posed (see analysis in Fox Reference Fox2019 and § 4) for material-density ratio
$Z = \rho _f / \rho _p \le 0.1$
.
In order to make the two-fluid model well posed for arbitrary
$Z$
, Fox, Laurent & Vié (Reference Fox, Laurent and Vié2020) included added mass as first proposed by Cook & Harlow (Reference Cook and Harlow1984). In the kinetic equation, the definition of a particle was extended to include some fluid (i.e. the added mass) travelling with the velocity of the particle. The effective material density of a particle
$\rho _e$
is defined using the sum of the particle mass and the added mass. As a result, the material-density ratio
$\rho _f / \rho _e$
remains finite when
$\rho _p \to 0$
. Besides added mass, Fox et al. (Reference Fox, Laurent and Vié2020) introduced an ad hoc model for the pfp-pressure tensor
${\unicode{x1D64B}}_{\!pfp}$
(or, equivalently, a model for
$- \boldsymbol {\Sigma }_{pfp}$
). They justified this additional hydrodynamic stress by the fact that, when
$Z \gg 1$
(e.g. bubbly flows), particle–particle interactions are strongly mediated by the fluid phase (Wang et al. Reference Wang, Yang, Zhang and Balachandar2021; Zhang Reference Zhang2021), as opposed to hard-sphere collisions. Nonetheless, while this augmented two-fluid model with added mass is well posed for arbitrary
$Z$
, the question of what is the exact form for
${\unicode{x1D64B}}_{\!pfp}$
still remains. Indeed, the general analysis by Zhang et al. (Reference Zhang, Ma and Rauenzahn2006) shows that the pfp-stress tensor should exist in any kinetic derivation of the two-fluid model. If this is true, then what is the exact expression for
${\unicode{x1D64B}}_{\!pfp}$
in the hard-sphere model of Fox (Reference Fox2019), and how is it related to
$\boldsymbol {F}_{\!\!fp}$
? Answering these questions is the principal objective of this work.
The remainder of this article is organised as follows. In § 2, we briefly review the BE kinetic model with emphasis on the particle-pair distribution function in the ideal-fluid–particle limit. The two-fluid model derived in Fox (Reference Fox2019) for an ideal fluid and spherical particles is also provided in § 2. The main result concerning the relationship between
$\boldsymbol {F}_{\!\!fp}$
and
${\unicode{x1D64B}}_{\!pfp}$
is derived in § 3. Finally, in § 4, we find explicit expressions for the eigenvalues of the two-fluid model for a constant-density fluid, and discuss how the results from this work can be applied to the two-fluid models containing additional physics discussed in Fox (Reference Fox2024).
2. Boltzmann–Enskog kinetic model
In Fox (Reference Fox2019), a hyperbolic two-fluid model was derived, starting from the BE kinetic equation, that includes the Archimedes force. The derivation used binary, hard-sphere collisions for particles with diameters
$d_1 \le d_2$
, and resulted in a two-fluid model for the limit
$d_1 \ll d_2$
, e.g. a molecular gas containing finite-size particles. In this section, we briefly review the principal assumptions and results for an ideal-fluid–particle flow.
2.1. Principal assumptions
For simplicity and clarity, we only consider the ideal-fluid case. The kinetic theory model in Fox (Reference Fox2019) uses the Enskog term for the two-particle distribution function at the contact point (i.e. on the surface of the particle with diameter
$d_p$
)
\begin{equation} f^{(2)}({\mathbf{v}}_f,{\mathbf{v}}_p) = f_f({\mathbf{v}}_f) f_p ({\mathbf{v}}_p) + \frac {1}{2} d_p f_p ({\mathbf{v}}_p) \, \mathbf{x}_{12} \boldsymbol {\cdot } \partial _{\mathbf{x}} f_f ({\mathbf{v}}_f) ,\end{equation}
where
$\mathbf{x}_{12}$
is the unit vector pointing from the surface contact point for collision with a fluid particle to the particle centre. The reader will recognise (2.1) as a Taylor-like expansion of the fluid-velocity distribution function (VDF), denoted by
$f_{\!f}$
, from the particle surface to its centre. The principal assumption used to derive the two-fluid model is that (2.1) is valid when spatial gradients are small enough to neglect higher-order terms.
Another assumption is that the VDF for an ideal fluid is Maxwellian
\begin{equation} f_{\!f} ({\mathbf{v}}_f) = \frac {\rho _f(t,\mathbf{x})}{[ 2 \pi \Theta _f (t,\mathbf{x}) ]^{3/2}} \exp \left ( - \frac { |{\mathbf{v}}_f - \boldsymbol{u}_f (t,\mathbf{x}) |^2}{2 \Theta _f (t,\mathbf{x})} \right ), \end{equation}
with velocity variance
$\Theta _f$
and fluid pressure
$p_f = \rho _f \Theta _f$
. Notice that
$f_{\!f}$
depends on
$\mathbf{x}$
and
$t$
through three moments, e.g. fluid density
$\rho _f$
, fluid velocity
$\boldsymbol{u}_f$
and fluid pressure
$p_f$
(one could also use
$\Theta _f$
in place of
$p_f$
, but
$p_f$
is more convenient for separating the Archimedes (buoyancy) force from the hydrodynamic force
$\boldsymbol {F}_{\!\!fp}$
.
$\Theta _f$
is proportional to the thermodynamic temperature
$T_f$
.) and not explicitly. Thus, the Enskog term in (2.1) can be written using the following expression for the spatial gradient (summation over
$i$
):
\begin{equation} \partial _{\mathbf{x}} f_{\!f} = \left (\frac {\partial f_{\!f} }{\partial \rho _f} \right ) \partial _{\mathbf{x}} \rho _f + \left (\frac {\partial f_{\!f} }{\partial u_{f,i}} \right ) \partial _{\mathbf{x}} u_{f,i} + \left (\frac {\partial f_{\!f} }{\partial p_f} \right ) \partial _{\mathbf{x}} p_f . \end{equation}
On the right-hand side of (2.3), only the partial derivatives of
$f_{\!f}$
depend on
${\mathbf{v}}_f$
through (2.2). In the kinetic theory model, the gradients of all hydrodynamic variables are considered, not just the pressure gradient needed for the Archimedes force.
Even before doing the detailed calculations to evaluate the collision integrals, we can already observe from (2.3) that the Archimedes force will come from the term involving
$\partial _{\mathbf{x}} p_f$
, while the two other terms generate
$\boldsymbol {F}_{\!pf}$
. In other words,
$\boldsymbol {F}_{\!pf}$
results from fluid density and velocity gradients across the particle radius. Note that the separation of
$f^{(2)}$
into two contributions is done in a more general formulation by Zhang et al. (Reference Zhang, Ma and Rauenzahn2006) and Zhang (Reference Zhang2021). The interest of using (2.1) as a specific example is that it allows us to find an exact expression for
$\boldsymbol {F}_{\!pf}$
, valid when spatial gradients of
$f_{\!f}$
are small. In contrast, Wang et al. (Reference Wang, Yang, Zhang and Balachandar2021) find an approximate correlation for
$\boldsymbol {F}_{\!pf}$
from particle-resolved DNS (direct-numerical simulation) that depends on particle-phase volume fraction
$\alpha _p$
and the particle Reynolds number
$Re_p$
. Generally speaking, small
$Re_p$
is a sufficient condition for small spatial variation of
$f_{\!f}$
. Thus, the main assumption from which all the results below follow is that (2.1) is a physically relevant approximation for
$f^{(2)}$
.
2.2. Collision integrals with Maxwellian VDF
The collision integral that generates the force of an ideal fluid on a spherical particle of diameter
$d_p$
uses the approximation in (2.1) for
$f^{(2)}$
. The first term on the right-hand side (
$f_{\!f} f_p$
) is the Boltzmann closure, while the second accounts for variations of
$f_{\!f}$
around the particle surface. The latter is the Enskog closure that keeps only the first derivative in
$\mathbf{x}$
. As shown in Fox (Reference Fox2019), the fluid-drag force depends on the first term in (2.1), while the Archimedes and hydrodynamic forces depend on the second. For clarity, we shall also assume that the particle VDF is Maxwellian
\begin{equation} f_p ({\mathbf{v}}_p) = \frac {\rho _p \alpha _p(t,\mathbf{x})}{[ 2 \pi \Theta _p (t,\mathbf{x}) ]^{3/2}} \exp \left ( - \frac {| {\mathbf{v}}_p - \boldsymbol{u}_p (t,\mathbf{x}) |^2}{2 \Theta _p (t,\mathbf{x})} \right ), \end{equation}
and the particle–particle collisions are elastic. A more general form where the particle velocity variance (
$\Theta _p$
is often referred to as granular temperature.)
$\Theta _p$
is a tensor (i.e. anisotropic Gaussian) and the collisions are inelastic can be found in Fox (Reference Fox2019). Note that
$\rho _p$
is constant and the fluid-phase volume fraction is defined by
$\alpha _f = 1 - \alpha _p$
.
The specific forms of the collision integral can be found in Fox (Reference Fox2019). For momentum, the forms (explicitly showing the dependence on
${\mathbf{v}}_f$
and
${\mathbf{v}}_p$
) are the vector
\begin{equation} C^{(0)} = \frac {3}{2 \rho _p d_p} \int _{\mathbb{R}^6} ({\mathbf{v}}_f - {\mathbf{v}}_p ) | {\mathbf{v}}_f - {\mathbf{v}}_p| f_{\!f}({\mathbf{v}}_f) f_p ({\mathbf{v}}_p) \, \textrm{d} {\mathbf{v}}_f \, \textrm{d} {\mathbf{v}}_p ,\end{equation}
and the vector
$C^{(1)}$
with components
\begin{equation} C^{(1)}_i = \frac {1}{5 \rho _p} \sum _j \int _{\mathbb{R}^6} I_{i,j} ({\mathbf{v}}_f - {\mathbf{v}}_p) f_p ({\mathbf{v}}_p) \partial _{x_j} f_{\!f} ({\mathbf{v}}_f) \, \textrm{d} {\mathbf{v}}_f \, \textrm{d} {\mathbf{v}}_p ,\end{equation}
where the symmetric second-order tensor
$I_{i,j}$
has diagonal elements
$I_{i,i} ({\mathbf{v}}) = {\mathbf{v}}^2 + 2 v_i^2$
and off-diagonal elements
$I_{i,j} ({\mathbf{v}}) = 2 v_i v_j$
(i.e.
$I_{i,j} ({\mathbf{v}}) = {\mathbf{v}}^2 \delta _{i,j} + 2 v_i v_j$
). This tensor is found by integration over the collision surface of a particle, and all such integrals have a closed form for integer velocity moments (see Fox Reference Fox2019 for details). Notice that
$I_{i,j}$
does not depend on
$\mathbf{x}$
or
$t$
, only on the independent variable
$\mathbf{v}$
. For a Maxwellian VDF, velocity moments of orders 0, 1 and 2 suffice to find the mass, momentum and total energy balances. The collision integrals for mass are null, and those for energy are given in Fox (Reference Fox2019). Here, we consider only (2.5) and (2.6) for momentum exchange.
2.2.1. Fluid-drag force
The parameter
$C^{(0)}$
represents the drag force and the integral is unclosed due to
$| {\mathbf{v}}_f - {\mathbf{v}}_p |$
. For the ideal-fluid case, both
$f_{\!f}$
and
$f_p$
are Maxwellian. Approximating the integrals yields an expression for the drag coefficient
$K$
. In general,
$\Theta _f \gg \Theta _p$
when
$Re_p$
is small so that
$K \propto 6 \, \alpha _p \rho _f \Theta _f^{1/2}/ d_p$
. This approximation is inaccurate and higher-order terms in the VDF are required (see van Beijeren & Dorfman Reference van Beijeren and Dorfman1980a
,Reference van Beijeren and Dorfman
b
for details). In practice,
$K$
is found using a drag correlation (Capecelatro Reference Capecelatro2022).
2.2.2. Forces due to spatial gradients
The parameter
$C^{(1)}$
represents all of the forces that depend on the spatial gradient
$\partial _{\mathbf{x}} f_{\!f}$
. Substituting (2.3) leads to three separate contributions (summation over
$j$
)
\begin{equation} \begin{gathered} \left [ \int _{\mathbb{R}^6} I_{i,j} ({\mathbf{v}}_f - {\mathbf{v}}_p) \left (\frac {\partial \ln f_{\!f} }{\partial \rho _f} \right ) f_p ({\mathbf{v}}_p) f_{\!f} ({\mathbf{v}}_f) \, \textrm{d} {\mathbf{v}}_f \, \textrm{d} {\mathbf{v}}_p \right ] \partial _{x_j} \rho _f \\ \left [ \int _{\mathbb{R}^6} I_{i,j} ({\mathbf{v}}_f - {\mathbf{v}}_p) \left (\frac {\partial \ln f_{\!f} }{\partial u_{f,k}} \right ) f_p ({\mathbf{v}}_p) f_{\!f} ({\mathbf{v}}_f) \, \textrm{d} {\mathbf{v}}_f \, \textrm{d} {\mathbf{v}}_p \right ] \partial _{x_j} u_{f,k} \\ \left [ \int _{\mathbb{R}^6} I_{i,j} ({\mathbf{v}}_f - {\mathbf{v}}_p) \left (\frac {\partial \ln f_{\!f} }{\partial p_f} \right ) f_p ({\mathbf{v}}_p) f_{\!f} ({\mathbf{v}}_f) \, \textrm{d} {\mathbf{v}}_f \, \textrm{d} {\mathbf{v}}_p \right ] \partial _{x_j} p_f ,\end{gathered} \end{equation}
where (the integral
$\int _{\mathbb{R}^3} |{\mathbf{v}}_f - \boldsymbol{u}_f |^2 f_{\!f} \, \textrm{d} {\mathbf{v}}_f$
(trace of velocity covariance) equals
$3 \rho _f \Theta _f = 3 p_f$
and
$\int _{\mathbb{R}^3} f_{\!f} \, \textrm{d} {\mathbf{v}}_f = \rho _f$
. Gaussian moments of order higher than two have known dependencies of
$\rho _f$
,
$\boldsymbol{u}_f$
and
$\Theta _f$
)
\begin{equation} \begin{gathered} \frac {\partial \ln f_{\!f} }{\partial \rho _f} = \frac {3}{2 \rho _f} - \frac {|{\mathbf{v}}_f - \boldsymbol{u}_f|^2}{6 p_f} , \ \frac {\partial \ln f_{\!f} }{\partial \boldsymbol{u}_f} = \frac {\rho _f ({\mathbf{v}}_f - \boldsymbol{u}_f)}{p_f} , \\ \frac {\partial \ln f_{\!f} }{\partial p_f} = - \frac {1}{2 p_f} + \frac {\rho _f |{\mathbf{v}}_f - \boldsymbol{u}_f|^2}{6 p_f^2} . \end{gathered} \end{equation}
An important technical point is that the three integrals in the square brackets in (2.7) are closed, i.e. they can be evaluated exactly and depend on
$\rho _f,\boldsymbol{u}_f,\Theta _f$
and
$\rho _p \alpha _p, \boldsymbol{u}_p, \Theta _p$
. (See Fox Reference Fox2019 for details.) Remarkably, the integral multiplying
$\partial _{x_j} p_f$
reduces to
$\alpha _p \delta _{i,j}$
, yielding the Archimedes force. The other two integrals yield (in Fox Reference Fox2019, an equivalent procedure is used to find the Archimedes and hydrodynamic forces. Because
$I_{i,j}$
does not depend on
$\mathbf{x}$
, the velocity integrals can be computed before taking the spatial derivatives of the fluid variables. This procedure reduces the number of integrals to evaluate from three to one, making it easier to extend to higher-order moments)
\begin{equation} \boldsymbol {F}_{\!pf} = \alpha _p ( \Theta _p {\unicode{x1D644}} + {\unicode{x1D64D}}) \boldsymbol {\cdot } \partial _{\mathbf{x}} \rho _f + \rho _f \alpha _p \frac {2}{5} [ \boldsymbol{u}_{fp} \boldsymbol {\cdot } (\partial _{\mathbf{x}} \boldsymbol{u}_{f}) + (\partial _{\mathbf{x}} \boldsymbol{u}_{f}) \boldsymbol {\cdot } \boldsymbol{u}_{fp} + ( \partial _{\mathbf{x}} \boldsymbol {\cdot } \boldsymbol{u}_f ) \boldsymbol{u}_{fp} ], \end{equation}
where
$\boldsymbol{u}_{fp} = \boldsymbol{u}_f - \boldsymbol{u}_p$
and
\begin{equation} {\unicode{x1D64D}} = \frac {1}{5} u_{fp}^2 {\unicode{x1D644}} + \frac {2}{5} \boldsymbol{u}_{fp} \otimes \boldsymbol{u}_{fp}, \end{equation}
with
$tr({\unicode{x1D64D}}) = u_{fp}^2$
. Note that, because here we separate out the term involving
$\Theta _p {\unicode{x1D644}}$
when defining
${\unicode{x1D64D}}$
, this definition is slightly different than in Fox (Reference Fox2019).
In Fox (Reference Fox2019), the exact expression for
$C^{(1)}$
for arbitrary VDFs (e.g. depending on the velocity covariance matrix, etc.) is computed. Thus, (2.9) represents the lowest-order approximation assuming Maxwellian VDFs. In any case, the accuracy of this result depends mainly on the accuracy of the model for
$f^{(2)}$
in (2.1).
2.3. Two-fluid model for fluid–particle flow
For completeness, we provide the two-fluid model with an ideal fluid in table 1. Our main objective in the next section is to demonstrate that the momentum-exchange term
$\boldsymbol {F}_{\!pf}$
in (2.9) provides a closed expression for the pfp-pressure tensor
${\unicode{x1D64B}}_{\!pfp}$
introduced in an ad hoc manner in Fox et al. (Reference Fox, Laurent and Vié2020). The tensor
${\unicode{x1D64D}}$
results from integrating over the collision surface of a spherical particle (see chapter 6 in Marchisio & Fox Reference Marchisio and Fox2013 for details). Notice that
${\unicode{x1D64D}}$
appears in the fluid-phase momentum flux (‘collisional’ flux) in the momentum balance in table 1 (see Fox Reference Fox2019 for details). This result is also found by volume averaging over the surface of spherical particles (Lhuillier Reference Lhuillier2023), and arises due to the excluded volume occupied by particles. This would suggest that the fluid-phase momentum flux in the two-fluid model with spherical particles should always depend on
${\unicode{x1D64D}}$
, regardless of how it is derived.
Table 1. Two-fluid model for flow of an ideal fluid and elastic hard-sphere particles with constant
$\rho _p$
.
${\boldsymbol {F}}_{\!pf} = - {\boldsymbol {F}}_{fp}$
is given in terms of the pfp-pressure tensor
${\unicode{x1D64B}}_{\!pfp}$
and the volume-average velocity
$\boldsymbol{u}_v = \alpha _p \boldsymbol{u}_p + \alpha _f \boldsymbol{u}_f$
in (3.5). To include added mass, as done in Fox et al. (Reference Fox, Laurent and Vié2020), it suffices to make the following substitutions:
$\alpha _p \to \alpha _p^\star$
,
$\rho _p \to \rho _e$
,
$\alpha _f \to \alpha _f^{\star}$
; and to include a mass balance for
$\alpha _p$
, as well as the added-mass exchange terms (see Boniou et al. Reference Boniou, Fox, Posey and Houim2024 for details). Note that, for clarity, the frictional pressure needed for dense granular flows (see, e.g. Houim & Oran Reference Houim and Oran2016) has not been included in the definition of
$p_p$
. With or without added mass, the frictional pressure depends only on
$\alpha _p$
.
3. Particle–fluid–particle pressure tensor
In this section, we show how to find a closed expression for the pfp-pressure tensor
${\unicode{x1D64B}}_{\!pfp}$
from the momentum-exchange term
$\boldsymbol {F}_{\!pf}$
in (2.9). The derivation of
${\unicode{x1D64B}}_{\!pfp}$
requires us to express the gradient of the fluid velocity
$\partial _{\mathbf{x}} \boldsymbol{u}_{f}$
in terms of the slip velocity
$\boldsymbol{u}_{pf}$
and the volume-average velocity
$\boldsymbol{u}_v$
defined by
\begin{equation} \boldsymbol{u}_v = \alpha _p \boldsymbol{u}_p + \alpha _f \boldsymbol{u}_f = \alpha _p \boldsymbol{u}_{pf} + \boldsymbol{u}_f. \end{equation}
From the mass balances in table 1, it can easily be shown that if the material densities (i.e.
$\rho _f$
and
$\rho _p$
) are constant, then
$\partial _{\mathbf{x}} \boldsymbol {\cdot } \boldsymbol{u}_v = 0$
. Thus,
$\boldsymbol{u}_v$
is the ‘natural’ choice for the reference velocity. The momentum-exchange terms in (2.9) can then be rewritten using the identity
\begin{equation} \partial _{\mathbf{x}} \boldsymbol{u}_{f} = \partial _{\mathbf{x}} \boldsymbol{u}_{v} + \partial _{\mathbf{x}} (\alpha _p \boldsymbol{u}_{fp}). \end{equation}
In particular, by rearranging the terms on the left-hand side below (which appears in
$\boldsymbol {F}_{\!\!fp}$
), we find that
\begin{align} \rho _f \alpha _p \frac {2}{5} [ \boldsymbol{u}_{fp} \boldsymbol {\cdot } (\partial _{\mathbf{x}} \boldsymbol{u}_{f}) + (\partial _{\mathbf{x}} \boldsymbol{u}_{f}) \boldsymbol {\cdot } \boldsymbol{u}_{fp} + ( \partial _{\mathbf{x}} \boldsymbol {\cdot } \boldsymbol{u}_f ) \boldsymbol{u}_{fp} ] \nonumber \\ = \rho _f \alpha _p \frac {2}{5} [ 2 \boldsymbol {\Gamma }_v + ( \partial _{\mathbf{x}} \boldsymbol {\cdot } \boldsymbol{u}_v ) {\unicode{x1D644}} ] \boldsymbol {\cdot } \boldsymbol{u}_{fp} + \rho _f \partial _{\mathbf{x}} \boldsymbol {\cdot } (\alpha _p^2 {\unicode{x1D64D}}), \end{align}
where
$\boldsymbol {\Gamma }_v = \frac {1}{2} [(\partial _{\mathbf{x}} \boldsymbol{u}_v) + (\partial _{\mathbf{x}} \boldsymbol{u}_v)^t]$
is the rate-of-deformation tensor for the volume-average velocity.
The pfp-pressure tensor is then defined by
\begin{equation} {\unicode{x1D64B}}_{\!pfp} = \rho _f \alpha _p^2 {\unicode{x1D64D}} ,\end{equation}
such that
$\partial _{\mathbf{x}} \boldsymbol {\cdot } {\unicode{x1D64B}}_{\!pfp} = \rho _f \partial _{\mathbf{x}} \boldsymbol {\cdot } (\alpha _p^2 {\unicode{x1D64D}}) + \alpha _p^2 {\unicode{x1D64D}} \boldsymbol {\cdot } \partial _{\mathbf{x}} \rho _f$
. The dependence of
${\unicode{x1D64B}}_{\!pfp}$
on
$\alpha _p^2$
is consistent with it arising from particle-pair interactions, while the dependence on
$\rho _f$
shows the role of the fluid. Using these results, the momentum-exchange term becomes
\begin{align} \boldsymbol {F}_{\!pf} = \partial _{\mathbf{x}} \boldsymbol {\cdot } {\unicode{x1D64B}}_{\!pfp} + \alpha _p \rho _f \frac {4}{5} \boldsymbol {\Gamma }_v \boldsymbol {\cdot } \boldsymbol{u}_{fp} + \alpha _p ( \Theta _p {\unicode{x1D644}} + \alpha _f {\unicode{x1D64D}}) \boldsymbol {\cdot } \partial _{\mathbf{x}} \rho _f + \alpha _p \rho _f \frac {2}{5} ( \partial _{\mathbf{x}} \boldsymbol {\cdot } \boldsymbol{u}_v ) \boldsymbol{u}_{fp} . \\[-16pt] \nonumber \end{align}
The first two terms on the right-hand side are present even when
$\rho _f$
is constant, in which case the last two are zero. The second term is a lift force due to gradients in
$\boldsymbol{u}_v$
(and not
$\boldsymbol{u}_f$
); however, in the dilute limit
$\alpha _p \to 0$
(i.e. the lift force on a single particle), these two velocities are nearly identical. Alternatively (i.e. using the traceless tensor
${\unicode{x1D64E}}_v$
as done in Fox Reference Fox2019), the tensor
$\boldsymbol {\Gamma }_{\!v}$
in the lift force can be made traceless and one third its trace combined with the final term in (3.5).
Given that
$\boldsymbol {F}_{\!\!fp} = - \boldsymbol {F}_{\!pf}$
, we can observe that by moving
$\partial _{\mathbf{x}} \boldsymbol {\cdot } {\unicode{x1D64B}}_{\!pfp}$
to the left-hand sides of the momentum balances in table 1,
${\unicode{x1D64B}}_{\!pfp}$
acts like a (positive) pressure term for the particle phase. For the fluid phase, we can combine the
${\unicode{x1D64B}}_{\!pfp}$
term with the momentum flux due to
${\unicode{x1D64D}}$
to find
\begin{equation} \partial _{\mathbf{x}} \boldsymbol {\cdot } [ \alpha _p \rho _f ( \Theta _p {\unicode{x1D644}} + {\unicode{x1D64D}} )] - \partial _{\mathbf{x}} \boldsymbol {\cdot } {\unicode{x1D64B}}_{\!pfp} = \partial _{\mathbf{x}} \boldsymbol {\cdot } [ \alpha _p \rho _f ( \Theta _p {\unicode{x1D644}} + \alpha _f {\unicode{x1D64D}} )], \end{equation}
and thus the combined term acts as a pressure as well. In agreement with Zhang (Reference Zhang2021), this would suggest that a pfp-pressure tensor should be present in every two-fluid model for disperse multiphase flows, regardless of how the model is derived.
4. Discussion and conclusions
4.1.
Physical interpretation of
$\boldsymbol {F}_{\!pf}$
The force
$\boldsymbol {F}_{\!pf}$
has two contributions: one due to the fluid-density gradient at constant pressure, and one due to the fluid-velocity gradient. The first is important in high-speed flows with shock waves. At a contact surface behind a shock wave, the fluid pressure is constant but the density gradient is large. Here,
$\boldsymbol {F}_{\!pf}$
provides an unsteady drag on the particles due to the density gradient (Boniou & Fox Reference Boniou and Fox2023). In low-speed flows, the fluid-density gradient is small (or zero), in which case
$\boldsymbol {F}_{\!pf}$
depends only on fluid-velocity gradients. Physically, when
$\boldsymbol{u}_f$
varies over the length scale
$d_p/2$
, a hydrodynamic force results. Because
$\boldsymbol{u}_f = \boldsymbol{u}_v + \alpha _p \boldsymbol{u}_{pf}$
and (with constant
$\rho _f$
)
$\partial _{\mathbf{x}} \boldsymbol {\cdot } \boldsymbol{u}_v =0$
, the hydrodynamic stress scales like
$\rho _f \alpha _p^2 u_{pf}^2$
(see Zhang et al. Reference Zhang, Ma and Rauenzahn2006; Zhang Reference Zhang2021 for details). Thus, in addition to the fluid-drag force found for uniform
$\boldsymbol{u}_f$
, there is ‘slip stress’ when
$\boldsymbol{u}_f$
varies across the particle radius (i.e. for non-zero
$\partial _{\mathbf{x}} \boldsymbol{u}_f$
). Finally, it is important to notice that
$\boldsymbol {F}_{\!pf}$
does not depend directly on gradients of
$\alpha _p$
. This is because at the scale of the particle,
$\alpha _p$
is not defined (or constant). Thus, in the kinetic model in (2.1) for an ideal-fluid–particle system, the Enskog contribution to
$f^{(2)}$
does not depend on the gradient of
$f_p$
. Nonetheless, recent work (Wang, Zhang & Balachandar Reference Wang, Zhang and Balachandar2024) has demonstrated that a volume-fraction gradient leads to a diffusion stress
$\boldsymbol {\Sigma }_d$
. The reader is referred to Wang et al. (Reference Wang, Zhang and Balachandar2024) for the exact definition, but it is important to note that such a stress will affect the hyperbolicity of the two-fluid model. Indeed, a similar term is widely used in two-fluid models for bubbly flows to make them well-posed (Lhuillier et al. Reference Lhuillier, Chang and Theofanous2013; Panicker et al. Reference Panicker, Passalacqua and Fox2018).
4.2. Is the two-fluid model well-posed?
As shown in Fox et al. (Reference Fox, Laurent and Vié2020), since the trace of
${\unicode{x1D64B}}_{\!pfp}$
is
$\rho _f \alpha _p^2 u_{pf}^2$
, the pfp-pressure will be large enough to keep the two-fluid model well posed. For an incompressible fluid (see table 2), the two eigenvalues for the particle phase
$\lambda _{1,2} = u_p + u_{fp} \hat {\lambda }_{1,2}$
can be computed analytically (Panicker et al. Reference Panicker, Passalacqua and Fox2018). Subtracting
$u_p=0$
and dividing the result by
$u_{fp}$
, we find the two scaled eigenvalues (this result is strongly dependent on using both the correct pfp-pressure and fluid-phase momentum flux. For example, using
$\partial _{\mathbf{x}} \boldsymbol {\cdot } {\unicode{x1D64B}}_{\!pfp}$
without
$\partial _{\mathbf{x}} \boldsymbol {\cdot } [ \alpha _p \rho _f ( \Theta _p {\unicode{x1D644}} + {\unicode{x1D64D}} )]$
yields very different (non-physical) eigenvalues for small
$\Theta _p$
)
\begin{equation} \hat {\lambda }_{1,2} = \frac {Z \alpha _p \alpha _f \pm X^{1/2}}{(\alpha _f + Z \alpha _p) \alpha _f} ,\end{equation}
with (if the particle–particle collisions are inelastic with coefficient of restitution
$e_c$
, then
$p_p = \alpha _p \Theta _p [1+2(1+e_c) \alpha _p g_0]$
. The second and third terms on the right-hand side of (4.2) then depend on
$e_c$
, but the qualitative behaviour of the eigenvalues remains the same)
\begin{equation} X = (1 - \hat {\Theta }_p) Z^2 \alpha _p^2 \alpha _f^2 + \frac {2}{\alpha _f} (3+\alpha _f) \alpha _p^2 Z \hat {\Theta }_p + ( \alpha _p^4 + 4 \alpha _f \alpha _p^2 +4 \alpha _p + 1) \hat {\Theta }_p ,\end{equation}
where
$\hat {\Theta }_p = \Theta _p / u_{fp}^2$
. As seen in figure 1, these eigenvalues are real and distinct for all
$\alpha _p$
when
$X \gt 0$
. As only the first term in (4.2) can be negative when
$\hat {\Theta }_p \gt 1$
, an ill-posed system may arise if
$Z$
is large enough. In the granular limit,
$Z=0$
and
$\hat {\lambda }_{1,2} \propto \pm \hat {\Theta }_p^{1/2}$
as expected, while for massless particles
$Z \to \infty$
and
$\hat {\lambda }_{1,2} = 1 \pm (1-\hat {\Theta }_p)^{1/2}$
. However, when added mass is included in the two-fluid model (see Fox et al. Reference Fox, Laurent and Vié2020; Boniou et al. Reference Boniou, Fox, Posey and Houim2024 for details), the effective density of the particles
$\rho _e$
depends on the added-mass coefficient
$C_m \approx 0.5$
(for the two-fluid model with added mass to be well posed for all
$\Theta _p$
, the minimum value is
$C_m \approx 0.1$
, which gives
$Z \approx 11$
when
$\rho _p \to 0$
) such that
$Z = \rho _f / \rho _e \le 3$
. In that case, the two eigenvalues are real for all
$\hat {\Theta }_p \ge 0$
. Likewise, when the frictional pressure is included in the particle-phase pressure
$p_p$
, the particle-phase eigenvalues have much larger magnitude when
$\alpha _p \gt 0.63$
(Houim & Oran Reference Houim and Oran2016). We can therefore conclude that the two-fluid model derived in Fox (Reference Fox2019) with the pfp-pressure tensor given in (3.4) and the added-mass model in Fox et al. (Reference Fox, Laurent and Vié2020), Boniou et al. (Reference Boniou, Fox, Posey and Houim2024) is well posed for all physically relevant conditions.
Table 2. One-dimensional two-fluid model with a constant-density fluid and constant granular temperature
$\Theta _p$
. The source terms are neglected because they do not affect the eigenvalues. The primitive variables are
$X = [ p_f, \alpha _p, u_p, u_f]$
where the fluid pressure is divided by
$\rho _f$
. Two eigenvalues for this system are
$\pm \infty$
(Panicker et al. Reference Panicker, Passalacqua and Fox2018), and the other two can be scaled to depend only on
$Z$
and
$\hat {\Theta }_p$
(Fox et al. Reference Fox, Laurent and Vié2020). In the particle-phase momentum balance, the second term is the granular pressure, the third term is the pfp contribution, while the fourth is due to the Archimedes force. The pfp and
$\alpha _p R$
contributions are combined in the fluid-phase momentum balance, yielding a positive ‘slip pressure’ as in (3.6).
Figure 1. Scaled particle-phase eigenvalues with an incompressible fluid found from (4.1) with (a)
$\hat {\Theta }_p =0$
and (b)
$\hat {\Theta }_p =1$
. Lines: solid black,
$Z = 1000$
; dotted green,
$Z=3$
; dash-dot blue,
$Z=1$
; dashed red,
$Z=0.001$
. With
$\hat {\Theta }_p =0$
, one scaled eigenvalue is zero. For larger
$\hat {\Theta }_p$
, the eigenvalues separate more rapidly with increasing
$\alpha _p$
. When added mass is included (Fox et al. Reference Fox, Laurent and Vié2020; Boniou et al. Reference Boniou, Fox, Posey and Houim2024), bubbly flow corresponds to
$Z \approx 3$
, and the eigenvalues are real for all
$\hat {\Theta }_p \ge 0$
.
Concerning the physical interpretation of the eigenvalues in figure 1, the results for
$\Theta _p = 0$
are of particular interest. In a frame of reference with
$u_p=0$
, the particles are fixed with zero mean and fluctuating velocity, and the slip velocity
$u_{fp}$
is in the positive direction (e.g. gravitational sedimentation). This results in
$\lambda _1 = 0$
so that no information about the particles travels upstream (i.e. ahead of the falling particles). On the other hand, the wave with
$\lambda _2 = 2 Z \alpha _p u_{fp} / (\alpha _f + Z \alpha _p)$
will propagate particle-phase mass and momentum downstream at speeds that can exceed the mean fluid velocity for
$Z \ge 1$
. This will cause a cloud of particles to spread only in the downstream (
$u_{fp} \gt 0$
) direction. In most applications (e.g. Boniou & Fox Reference Boniou and Fox2023),
$0 \lt \hat {\Theta }_p \ll 1$
so that information propagates in both directions. Nonetheless, due to its dependence on
$Z$
, including the pfp-pressure term (and added mass) is crucial for lightweight particles such as bubbles (Risso Reference Risso2018) or even inertial granular suspensions (Guazzelli & Pouliquen Reference Guazzelli and Pouliquen2018).
In summary, like the Euler equations for ideal-fluid flow, the two-fluid model in table 1 for ideal-fluid–particle flow is the simplest possible model containing the minimum physics needed to be well posed. Thus, it can be used as the starting point for adding more physics such as viscous effects, mass transfer and chemical reacting species (Boniou et al. Reference Boniou, Fox, Posey and Houim2024). Unlike the ‘standard’ two-fluid model used in most academic and commercial codes, the model in table 1 has the fluid-phase flux term
${\unicode{x1D64D}}$
and the pfp-pressure tensor
${\unicode{x1D64B}}_{\!pfp}$
, both of which arise (like the Archimedes force) from the excluded volume of the particles. Given this fact, it is reasonable to expect that the fluid–particle model in table 1 can be used as a starting point for other particle shapes besides spheres, albeit with corrections to some of the exchange terms (e.g. the drag coefficient
$K$
).
Acknowledgements.
The author would like to thank the Massot group at CMAP–École Polytechnique for motivating him to take a fresh look at the relationship between
$\boldsymbol {F}_{\!pf}$
and
${\unicode{x1D64B}}_{\!pfp}$
. The insightful comments and comprehensive feedback from the anonymous referees was highly appreciated by the author and greatly improved the readability and focus of this work.
Funding.
For the performance of this research, the author was supported by the Fulbright–Tocqueville Distinguished Scholar and Jean D’Alembert Senior Professor Chair at the Université de Paris-Saclay.
Declaration of interests.
The author reports no conflict of interest.
Open access
