2.1. Derivation of two-fluid momentum balances and definition of buoyancy closure

In this subsection, we summarize the results of the mathematically rigorous derivation (i.e. no approximations involved) by Pähtz et al. (Reference Pähtz, Chen, Zhu, Tholen and He2026) of the macroscopic two-fluid momentum balances from averaging the microscopic equations of motion of the fluid and nearly rigid particles. However, note that their assumption of nearly rigid particles (clarified shortly) is not required for the buoyancy-related analysis we perform afterwards, since considering deformable and/or compressible particles, such as droplets and bubbles, has no effect on the average fluid phase momentum balance and the micromechanical expressions for the quantities therein as long as interfacial mass transfer does not occur (Fintzi & Pierson Reference Fintzi and Pierson2026).

Pähtz et al. (Reference Pähtz, Chen, Zhu, Tholen and He2026) defined the microscopic system as follows: the fluid–particle mixture is contained within the domain $\mathbb{V}_\infty =\mathbb{V}_{\mathit{f}}[t]\cup \mathbb{V}_{\mathit{s}}[t]$ composed of the fluid phase $\mathbb{V}_{\mathit{f}}$ and dispersed phase $\mathbb{V}_{\mathit{s}}$ . The latter is subdivided into individual particles $p$ : $\mathbb{V}_{\mathit{s}}[t]=\cup _p\mathbb{V}^p[t]$ . The corresponding indicator functions $X_\infty [\boldsymbol{x}]$ , $X_{\mathit{f}}[\boldsymbol{x},t]$ , $X_{\mathit{s}}[\boldsymbol{x},t]$ and $X^p[\boldsymbol{x},t]$ , with $\boldsymbol{x}\in \mathbb{R}^3$ the spatial coordinate, are defined as being equal to $1$ within the interiors and equal to $0$ outside of $\mathbb{V}_\infty$ , $\mathbb{V}_{\mathit{f}}$ , $\mathbb{V}_{\mathit{s}}$ and $\mathbb{V}^p$ , respectively. The fluid satisfies general microscopic mass and momentum balance equations,

(2.1) \begin{align} \partial _t\rho +\boldsymbol{\nabla }\boldsymbol{\cdot }\rho \boldsymbol{u}&=0,\\[-12pt]\nonumber \end{align}

(2.2) \begin{align} \rho {\mathrm{D}}_t\boldsymbol{u}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }+\rho \boldsymbol{g}, \end{align}

where ${\mathrm{D}}_t\equiv \partial _t+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }$ . In contrast to before, $\boldsymbol{\sigma }$ is a stress tensor associated with an arbitrary fluid flow rheology, and the density $\rho$ is not necessarily constant and therefore distinguished from the previous $\rho _{\mathit{f}}$ . In particular, the microscopic fields $\rho [\boldsymbol{x},t]$ and $\boldsymbol{u}[\boldsymbol{x},t]$ are not only defined for $\boldsymbol{x}\in \mathbb{V}_{\mathit{f}}$ but also for $\boldsymbol{x}\in \mathbb{V}^p$ , where they coincide with a particle $p$ ’s (not necessarily constant) interior mass density distribution and its velocity $\boldsymbol{v}^p_{\uparrow} +\boldsymbol{\omega }^p\times \boldsymbol{r}^p$ , respectively. Here $\boldsymbol{v}^p_{\uparrow} [t]$ is its centre-of-mass velocity, $\boldsymbol{\omega }^p[t]$ its angular velocity and $\boldsymbol{r}^p[\boldsymbol{x},t]\equiv \boldsymbol{x}-\boldsymbol{x}^p[t]$ the coordinate relative to its centre-of-mass location $\boldsymbol{x}^p[t]$ . A density field $\rho _\ast$ based on the average mass density $M^p/V^p$ of each particle $p$ , with $M^p$ its mass and $V^p$ its volume, and a velocity field $\boldsymbol{u}_{\uparrow}$ without particle-rotational contributions are also defined as $\rho _\ast \equiv X_{\mathit{f}}\rho +\sum _pX^pM^p/V^p$ and $\boldsymbol{u}_{\uparrow} \equiv X_{\mathit{f}}\boldsymbol{u}+\sum _pX^p\boldsymbol{v}^p_{\uparrow}$ , respectively. The particles are non-Brownian and no interfacial mass transfer is allowed to occur (impermeable boundary conditions), implying that the normal velocity $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}^p$ is continuous at the surface $\mathbb{S}^p=\partial \mathbb{V}^p$ , where $\boldsymbol{n}^p$ is the corresponding unit normal vector pointing outward the particle $p$ . They experience the instantaneous, local hydrodynamic force $\boldsymbol{F}^p_{\mathit{h}}=\int _{\mathbb{S}^p}\boldsymbol{n}^p\boldsymbol{\cdot }\boldsymbol{\sigma }{\mathrm{d}} S$ , which defines the microscopic hydrodynamic force density field

(2.3) \begin{align} \boldsymbol{f}_{\!\mathit{h}}\equiv \sum \nolimits _p\frac {X^p}{V^p}\boldsymbol{F}^p_{\mathit{h}}=\sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{S}^p}\boldsymbol{n}^p\boldsymbol{\cdot }\boldsymbol{\sigma }{\mathrm{d}} S. \end{align}

Furthermore, they interact via contacts, where slight transient elastic deformations are permitted (i.e. particles are not fully rigid). However, it is assumed that these do not alter the particles’ overall shapes significantly and that the contact area of a contacting particle pair $pq$ can be represented by a single contact point $\boldsymbol{x}^{pq}_{\mathit{c}}=\boldsymbol{x}^{qp}_{\mathit{c}}$ . The total contact force acting on a particle $p$ is denoted as $\boldsymbol{F}^p_{\mathit{c}}=\sum _q\boldsymbol{F}^{pq}_{\mathit{c}}$ , resulting from interactions $\boldsymbol{F}^{pq}_{\mathit{c}}=-\boldsymbol{F}^{qp}_{\mathit{c}}$ with other particles $q$ (by convention, $\boldsymbol{F}^{pp}_{\mathit{c}}=0$ ).

The fluid’s and particles’ microscopic equations of motion are averaged using an averaging procedure $\langle \boldsymbol{\cdot }\rangle$ that satisfies the following rules (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026):

(2.4) \begin{align} &\langle cA_1+A_2\rangle =c\langle A_1\rangle +\langle A_2\rangle ,&&\langle \boldsymbol{\nabla }A\rangle =\boldsymbol{\nabla }\langle A\rangle ,&&\langle \partial _tA\rangle =\partial _t\langle A\rangle ,&&\langle AX_\infty \rangle =\langle A\rangle , \end{align}

where the $A$ -s are arbitrary fields well-defined for $\boldsymbol{x}\in \mathbb{V}_\infty$ and $c$ is a constant. In addition, the Reynolds rule,

(2.5) \begin{align} \langle \langle A_1\rangle A_2\rangle =\langle A_1\rangle \langle A_2\rangle , \end{align}

must also be obeyed. It ensures uniqueness of the averaged fields and, when excluding the trivial (and generally unphysical) cases of the identity operator and operators that average over all times and/or space, essentially states that the averaging should include the entire ensemble of microscopic (i.e. fluctuating) states contributing to a macroscopic (i.e. average) state at a given $\boldsymbol{x}$ and $t$ , while not intermingling different macroscopic states (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026). At a mathematically exact level, this rule is generally violated for spatiotemporal averaging procedures, since they intermingle generally different macroscopic states from space and time surrounding $\boldsymbol{x}$ and $t$ . Consider, for example, a vertically stratified fluid at rest whose properties depend only on the vertical coordinate $z$ , where all microscopic states at a given $z$ are identical to each other and the corresponding macroscopic state. Then, spatial averaging that involves the $z$ -direction will intermingle different macroscopic states and cause differences between averaged and non-averaged quantities, even though physical intuition dictates that such differences should not occur for a system at rest.

The only spatiotemporal averaging procedures that identically satisfy all the rules in (2.4) and (2.5) are the infinite time average in statistically steady systems and infinite averages over spatial coordinates in which the system is homogeneous. If the ergodic hypothesis is true, they are equivalent to ensemble averaging, which also satisfies these rules identically.

Based on $\langle \boldsymbol{\cdot }\rangle$ , phase averages and mass-weighted phase averages are also defined (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026),

(2.6) \begin{align} \langle A\rangle ^{\,\mathit{f}}&\equiv \langle X_{\mathit{f}}A\rangle /\beta _{\mathit{f}},&\langle A\rangle ^{\mathit{s}}&\equiv \langle X_{\mathit{s}}A\rangle /\beta _{\mathit{s}},&\langle A\rangle &=\beta _{\mathit{f}}\langle A\rangle ^{\,\mathit{f}}+\beta _{\mathit{s}}\langle A\rangle ^{\mathit{s}},\\[-12pt]\nonumber \end{align}

(2.7) \begin{align} \langle A\rangle ^{\,\mathit{f}}_{\rho _\ast }&\equiv \langle \rho _\ast A\rangle ^{\,\mathit{f}}/\rho _{\mathit{f}},&\langle A\rangle ^{\mathit{s}}_{\rho _\ast }&\equiv \langle \rho _\ast A\rangle ^{\mathit{s}}/\rho _{\mathit{s}},&\langle A\rangle _{\rho _\ast }&\equiv \rho _{\mathit{m}}^{-1}\big(\beta _{\mathit{f}}\rho _{\mathit{f}}\langle A\rangle ^{\,\mathit{f}}+\beta _{\mathit{s}}\rho _{\mathit{s}}\langle A\rangle ^{\mathit{s}}\big), \end{align}

where $\beta _{\mathit{f}}\equiv \langle X_{\mathit{f}}\rangle$ is the fluid phase volume fraction, $\beta _{\mathit{s}}\equiv \langle X_{\mathit{s}}\rangle$ the dispersed-phase volume fraction, $\rho _{\mathit{f}}\equiv \langle \rho _\ast \rangle ^{\,\mathit{f}}$ the average fluid density, $\rho _{\mathit{s}}\equiv \langle \rho _\ast \rangle ^{\mathit{s}}$ the average dispersed-phase density, and $\rho _{\mathit{m}}\equiv \langle \rho _\ast \rangle =\beta _{\mathit{f}}\rho _{\mathit{f}}+\beta _{\mathit{s}}\rho _{\mathit{s}}$ the average mixture density. Note that these definitions are consistent with our prior usage of these quantities. In particular, the non-averaged value of $\rho _{\mathit{f}}$ was considered constant and was therefore identical to the averaged value defined here. Furthermore, note that the use of dispersed-phase averages for nearly rigid particles, which are weighted by the indicator function $X^p$ rather than the conventional multivariate delta distribution $\delta [\boldsymbol{x}-\boldsymbol{x}^p]$ , generally lead to particle-rotational terms in translational dispersed-phase balance equations (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026), see below.

With the above definitions of the microscopic system and averaged quantities, the mathematically exact averaged two-fluid momentum balances can be written as (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026)

(2.8) \begin{align} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}},\\[-12pt]\nonumber \end{align}

(2.9) \begin{align} \beta _{\mathit{s}}\rho _{\mathit{s}}{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}+\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{g}+\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}}, \end{align}

where $\boldsymbol{u}_{\mathit{f}}\equiv \langle \boldsymbol{u}_{\uparrow} \rangle _{\rho _\ast }^{\,\mathit{f}}$ is the average fluid phase velocity, $\boldsymbol{u}_{\mathit{s}}\equiv \langle \boldsymbol{u}_{\uparrow} \rangle _{\rho _\ast }^{\mathit{s}}$ the average dispersed-phase velocity without particle-rotational contributions, ${\mathrm{D}}^{\,\mathit{f}}_t\equiv \partial _t+\boldsymbol{u}_{\mathit{f}}\boldsymbol{\cdot }\boldsymbol{\nabla }$ and ${\mathrm{D}}^{\mathit{s}}_t\equiv \partial _t+\langle \boldsymbol{u}\rangle _{\rho _\ast }^{\mathit{s}}\boldsymbol{\cdot }\boldsymbol{\nabla }$ are the associated material derivatives, and $\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}\equiv -\beta _{\mathit{f}}\rho _{\mathit{f}}\langle (\boldsymbol{u}_{\uparrow} -\boldsymbol{u}_{\mathit{f}} ) (\boldsymbol{u}_{\uparrow} -\boldsymbol{u}_{\mathit{f}} )\rangle _{\rho _\ast }^{\,\mathit{f}}$ and $\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}\equiv -\beta _{\mathit{s}}\rho _{\mathit{s}}\langle (\boldsymbol{u}-\langle \boldsymbol{u}\rangle _{\rho _\ast }^{\mathit{s}} ) (\boldsymbol{u}_{\uparrow} -\boldsymbol{u}_{\mathit{s}} )\rangle _{\rho _\ast }^{\mathit{s}}$ the associated Reynolds stresses, respectively. Note that, when all particles are sphere-symmetrical, $\langle \boldsymbol{u}\rangle _{\rho _\ast }^{\mathit{s}}$ can be formally replaced by $\boldsymbol{u}_{\mathit{s}}$ within the definitions of ${\mathrm{D}}^{\mathit{s}}_t$ and $\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}$ , since their particle-rotational contributions cancel each other out. Furthermore, $\boldsymbol{\sigma ^{\,\mathit{f}}}\equiv \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ and $\boldsymbol{\sigma ^{\mathit{s}}}$ are the effective mechanical fluid phase and dispersed-phase stress tensors, respectively. (Note that this notation differs from Pähtz et al. (Reference Pähtz, Chen, Zhu, Tholen and He2026), who included the Reynolds stresses $\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}$ and $\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}$ within the definitions of $\boldsymbol{\sigma ^{\,\mathit{f}}}$ and $\boldsymbol{\sigma ^{\mathit{s}}}$ , respectively.) The detailed micromechanical expressions for $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ and $\boldsymbol{\sigma ^{\mathit{s}}}$ are provided in Appendix A, but they are not relevant for the coming analysis. The only relevant aspects are that the pseudostress $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ emerges due to the traction by the stress $\boldsymbol{\sigma }$ on the particles’ surfaces and that it satisfies (Supplementary material of Pähtz et al. (Reference Pähtz, Chen, Zhu, Tholen and He2026))

(2.10) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}=\beta _{\mathit{s}}\left \langle \sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{S}^p}\boldsymbol{n}^p\boldsymbol{\cdot }\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}{\mathrm{d}} S\right \rangle ^{\mathit{s}}+\langle \boldsymbol{\sigma }^{\textrm{T}}\rangle ^{\,\mathit{f}}\boldsymbol{\cdot }\boldsymbol{\nabla }\beta _{\mathit{s}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}, \end{align}

where $\boldsymbol{\cdot }^{\textrm{T}}$ denotes the transpose and $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}$ is the pseudostress that emerges due to the traction by the stress $\boldsymbol{\sigma }-\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ on the particles’ surfaces (Appendix A). The latter pseudostress therefore vanishes for fluid–particle systems at complete rest due to $\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}=\boldsymbol{\sigma }$ .

The averaged two-fluid momentum balances (2.8) and (2.9) contain various terms that require closure. However, throughout most of this paper (except in § 5.1), we do not employ any closure but the buoyancy closure, defined as the expression for $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ that separates generalized-drag contributions $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ from the total interfacial force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}}$ :

(2.11) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}}=\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}+\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}. \end{align}

The goal is to do this in a manner that is mathematically consistent with the exact micromechanical expressions for the other terms in (2.8) and does not cause disagreements with the particle-resolving simulations and thought experiment presented in § 3. In particular, to be consistent with the micromechanical definition of $\boldsymbol{f}_{\!\mathit{h}}$ , (2.3), it should be possible to express $\boldsymbol{f}_{\!\mathit{B}}$ in an analogous manner as

(2.12) \begin{align} \boldsymbol{f}_{\!\mathit{B}}=\sum \nolimits _p\frac {X^p}{V^p}\left \langle \boldsymbol{F}^p_{\mathit{B}}\right \rangle =\sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{S}^p}\boldsymbol{n}^p\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }{\mathrm{d}} S, \end{align}

where the average generalized-buoyancy force on a particle $p$ , $\langle \boldsymbol{F}^p_{\mathit{B}}\rangle$ , results from the integrated traction by some kind of macroscopic stress $\boldsymbol{\sigma _\ast }$ (i.e. a stress that satisfies $\langle \boldsymbol{\sigma _\ast }\rangle =\boldsymbol{\sigma _\ast }$ ) over its surface $\mathbb{S}^p$ . In two-fluid continuum models, one usually assumes that macroscopic stress divergences such as $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }$ do not vary much within a particle $p$ ’s volume $\mathbb{V}^p$ (Jackson Reference Jackson2000), in which case $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ approximates as

(2.13) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}&=\beta _{\mathit{s}}\left \langle \sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{S}^p}\boldsymbol{n}^p\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }{\mathrm{d}} S\right \rangle ^{\mathit{s}}=\left \langle X_{\mathit{s}}\sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{V}^p}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }{\mathrm{d}} V\right \rangle \nonumber\\[5pt] &=\left \langle \sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{V}^p}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }{\mathrm{d}} V\right \rangle \approx \left \langle \sum \nolimits _pX^p\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }\right \rangle =\beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }, \end{align}

where we used $X_{\mathit{s}}X^p=X^p$ and (2.6). We call this the ‘small-particle approximation’ of a buoyancy closure. However, buoyancy closures in the small-particle approximation lead to an inconsistency for fluid–particle systems at complete rest, as shown in the next subsection.

2.2. Thought experiment: fluid–particle systems at complete rest

For a fluid–particle system at complete rest, $\rho _{\mathit{f}}=\rho$ and $\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}=\boldsymbol{\sigma }$ at locations within the fluid phase, $\boldsymbol{x}\in V_{\mathit{f}}$ , and therefore $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}=0$ . Likewise, velocities and generalized-drag forces vanish due to the absence of motion. Then, the generalized-buoyancy force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ must be equal to the hydrodynamic force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}}$ , which implies $\boldsymbol{\sigma _\ast }=\boldsymbol{\sigma }$ for $\boldsymbol{x}\in V_{\mathit{f}}$ from comparing (2.12) with (2.3). Using $\boldsymbol{\sigma ^{\,\mathit{f}}}\equiv \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ and (2.10), it follows that the macroscopic fluid phase momentum balance (2.8) reduces to the stationary microscopic fluid phase momentum balance,

(2.14) \begin{align} X_{\mathit{f}}(\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }+\rho \boldsymbol{g})=0, \end{align}

as required. However, this consistency property is lost when using a buoyancy closure in the small-particle approximation. Then, (2.14) is recovered generally only for fluids where $\rho$ does not vary with $\boldsymbol{x}\in V_{\mathit{f}}$ . In particular, a particle bed immersed in a density-stratified fluid at rest would violate this property. This is not just a minor nuisance but is indirect evidence of a more fundamental issue. As we show in § 5.1, it is precisely the small-particle approximation that makes buoyancy cause ill-posedness of some two-fluid continuum models.

2.3. Previously proposed buoyancy closures

When substituting (2.10) via $\boldsymbol{\sigma ^{\,\mathit{f}}}\equiv \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ in the fluid phase momentum balance (2.8) and using the definition of the hydrodynamic microscopic force density $\boldsymbol{f}_{\!\mathit{h}}$ , (2.3), one obtains

(2.15) \begin{align} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}=\beta _{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}\big)+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\!\left \langle \sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{S}^p}\!\boldsymbol{n}^p\boldsymbol{\cdot }(\boldsymbol{\sigma }-\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}){\mathrm{d}} S\right \rangle ^{\mathit{s}}\!. \end{align}

This mathematical structure led Zhang et al. (Reference Zhang, VanderHeyden, Zou and Padial-Collins2007) to indirectly, and later others (e.g. Wang et al. Reference Wang, Zhang and Balachandar2024) to directly, suggest that the last term on the right-hand side of (2.15) is the negative generalized-drag force density, $-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ , and that

(2.16)

is therefore the generalized-buoyancy force density. This closure expression implies that $\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ is the background flow stress responsible for buoyancy, that is, $\boldsymbol{\sigma _\ast }=\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ in (2.12).

Ishii (Reference Ishii1975) proposed an alternative closure. They started with a different form of the fluid phase momentum balance,

(2.17) \begin{align} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}=\beta _{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}+\big\langle \big(\boldsymbol{\sigma }^{\textrm{T}}-\langle \boldsymbol{\sigma }^{\textrm{T}}\rangle ^{\,\mathit{f}}\big)\boldsymbol{\cdot }\boldsymbol{\nabla }X_{\mathit{s}}\big\rangle , \end{align}

which follows from (2.8) and the identity $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}}=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }^{\,\mathit{f}}_{\mathit{s}}-\left \langle \boldsymbol{\sigma }^{\textrm{T}}\boldsymbol{\cdot }\boldsymbol{\nabla }X_{\mathit{s}}\right \rangle$ (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026). They then interpreted the last term on the right-hand side of (2.17) as the negative generalized-drag force density, $-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}=\langle (\boldsymbol{\sigma }^{\textrm{T}}-\langle \boldsymbol{\sigma }^{\textrm{T}}\rangle ^{\,\mathit{f}})\boldsymbol{\cdot }\boldsymbol{\nabla }X_{\mathit{s}}\rangle$ . From comparison with (2.15), one can see that this interpretation is equivalent to assuming

(2.18)

However, this closure does not allow the identification of a background flow stress $\boldsymbol{\sigma _\ast }$ that is responsible for buoyancy. This is not surprising, since Ishii (Reference Ishii1975), as it was the fashion at the time (e.g. Drew Reference Drew1983), falsely treated the terms $\langle \boldsymbol{\sigma }^{\textrm{T}}\boldsymbol{\cdot }\boldsymbol{\nabla }X_{\mathit{s}}\rangle$ , $\langle \langle \boldsymbol{\sigma }^{\textrm{T}}\rangle ^{\,\mathit{f}}\boldsymbol{\cdot }\boldsymbol{\nabla }X_{\mathit{s}}\rangle$ , and $\langle (\boldsymbol{\sigma }^{\textrm{T}}-\langle \boldsymbol{\sigma }^{\textrm{T}}\rangle ^{\,\mathit{f}})\boldsymbol{\cdot }\boldsymbol{\nabla }X_{\mathit{s}}\rangle$ as actual fluid–particle interaction force densities, even though their micromechanical definitions do not take the particles’ surface geometries into account. In addition, even in the small-particle approximation, in which case (2.18) becomes $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}\approx \beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}$ , it is still not possible to identify a background flow stress $\boldsymbol{\sigma _\ast }$ , since this would require an expression of the form $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}\approx \beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma _\ast }$ as shown in (2.13). In other words, the closure (2.18) fails at a conceptual level.

Two further buoyancy closures based on the small-particle approximation have been put forward. Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013) essentially proposed

(2.19)

since they assumed that a mixture-viscosity-based shear stress, which micromechanically originates from the stress $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}$ (Jackson Reference Jackson1997), contributes to the background flow stress responsible for buoyancy. Furthermore, Anderson & Jackson (Reference Anderson and Jackson1967) and Jackson (Reference Jackson2000), though their micromechanical expressions for the stress tensor $\boldsymbol{\sigma ^{\,\mathit{f}}}$ were based on different levels of approximations, proposed the same buoyancy closure, henceforth called the ‘Jackson closure’:

(2.20) \begin{align} \text{Jackson closure}:\qquad \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}=\beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\big(\boldsymbol{\sigma ^{\,\mathit{f}}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}\big). \end{align}

Noting that $\boldsymbol{\sigma ^{\,\mathit{f}}}\equiv \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}\approx \langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}$ in the small-particle approximation, (2.19) and (2.20) differ from $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}\approx \beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ , the small-particle approximation of (2.16), in the additional inclusion of $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}$ and $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}$ , respectively, within the divergence operator. We comment on these closures in the next subsection.

2.4. Physically inspired buoyancy criterion

We propose that a physically meaningful buoyancy closure should satisfy the criterion

(2.21) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}=\beta _{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}, \end{align}

with a generalized-buoyancy force density that respects the geometry of the particles and therefore the mathematical structure (2.12)

(2.22) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}=\beta _{\mathit{s}}\left \langle \sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{S}^p}\boldsymbol{n}^p\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}{\mathrm{d}} S\right \rangle ^{\mathit{s}}. \end{align}

When integrating (2.21) over an infinitesimally small mixture volume ${\mathrm{d}} V$ , it states that the net flux ${\mathrm{d}} V\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}$ of average fluid momentum into ${\mathrm{d}} V$ decomposes into the net flux ${\mathrm{d}} V_{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ of average fluid momentum into ${\mathrm{d}} V_{\mathit{f}}\equiv \beta _{\mathit{f}}{\mathrm{d}} V$ , the average fluid phase portion of ${\mathrm{d}} V$ and the net flux ${\mathrm{d}} V_{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ of average fluid momentum into ${\mathrm{d}} V_{\mathit{s}}\equiv \beta _{\mathit{s}}{\mathrm{d}} V$ , the average dispersed-phase portion of ${\mathrm{d}} V$ . Importantly, both momentum flux contributions must be caused by the same stress $\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ , which can then be interpreted as the average background flow stress.

Using $\beta _{\mathit{f}}\langle \rho {\mathrm{D}}_t\boldsymbol{u}\rangle ^{\,\mathit{f}}=\beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}$ (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026) and (2.11), the criterion (2.21) allows transforming the fluid phase momentum balance (2.8) into the form

(2.23) \begin{align} \beta _{\mathit{f}}\langle \rho {\mathrm{D}}_t\boldsymbol{u}\rangle ^{\,\mathit{f}}=\beta _{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}. \end{align}

Its mathematical structure differs from the fluid phase momentum balance resulting from the closure by Zhang et al. (Reference Zhang, VanderHeyden, Zou and Padial-Collins2007), (2.15), in a single aspect: (2.23) does not contain a structural analogue to the term $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}$ in (2.15) on its right-hand side. When integrated over an infinitesimally small mixture volume ${\mathrm{d}} V$ , (2.23) states that the average temporal change ${\mathrm{d}} V_{\mathit{f}}\langle \rho {\mathrm{D}}_t\boldsymbol{u}\rangle ^{\,\mathit{f}}$ of fluid momentum within the pure-fluid volume ${\mathrm{d}} V_{\mathit{f}}$ is due to the net flux ${\mathrm{d}} V_{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ of average fluid momentum into ${\mathrm{d}} V_{\mathit{f}}$ and the body force, consisting of the gravitational and generalized-drag force terms. This is completely analogous to the microscopic fluid phase momentum balance (2.2) integrated over an infinitesimally small pure-fluid volume ${\mathrm{d}} V_{\mathit{f}}$ , where the temporal change ${\mathrm{d}} V_{\mathit{f}}\rho {\mathrm{D}}_t\boldsymbol{u}$ of fluid momentum within ${\mathrm{d}} V_{\mathit{f}}$ is due to the net flux ${\mathrm{d}} V_{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma }$ of fluid momentum into ${\mathrm{d}} V_{\mathit{f}}$ and the body force, consisting of only the gravitational term in this case. By contrast, such an analogy would not exist if there were mechanical stress divergence terms in the fluid phase momentum balance that were not being multiplied by $\beta _{\mathit{f}}$ , such as the term $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}$ in (2.15). Such terms would then still encode net fluxes of average fluid momentum into the mixture, as opposed to only the average fluid phase portion of the mixture, and therefore into the particles. It would mean that the net momentum transfer from the fluid to the particles would not be completely ‘outsourced’ to the interfacial force density terms, contrary to what we propose a buoyancy closure should actually achieve. In fact, in (2.23), fluid-to-particle momentum transfer is entirely encoded in the generalized-drag force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ .

The key questions now are: Does a stress $\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ that satisfies the criterion (2.21) exist? And, if it exists, what is its micromechanical expression? These questions will be positively answered in § 4. At the present stage, one can already infer from comparing (2.23) with (2.15) that $\boldsymbol{\hat \sigma ^{\,\mathit{f}}}=\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ for fluid–particle systems at complete rest due to $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}=0$ . By contrast, generally $\boldsymbol{\sigma ^{\,\mathit{f}}}\ne \langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ for such systems due to (2.10). Furthermore, from (2.21), it can already be seen that $\boldsymbol{\hat \sigma ^{\,\mathit{f}}}\approx \boldsymbol{\sigma ^{\,\mathit{f}}}$ when applying the small-particle approximation $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}\approx \beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ . As a consequence, (2.22) reduces to the closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), in this approximation.

Finally, it is worth noting that, when altering the criterion (2.21) in a manner that leads to an alternative version of (2.23) in which the average temporal change of fluid momentum per unit volume, $\beta _{\mathit{f}}\langle \rho {\mathrm{D}}_t\boldsymbol{u}\rangle ^{\,\mathit{f}}$ , on its left-hand side is replaced by the temporal change of the average fluid momentum per unit volume, $\beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}$ , one would obtain the Jackson closure, (2.20), in the small-particle approximation, which additionally includes the Reynolds stress $\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}$ within the divergence operator. However, this appears inconsistent with the micromechanical expressions for $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}}$ and $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ , based on (2.3) and (2.12), respectively, since they refer to the average force per unit volume acting on actual particles rather than the force per unit volume acting on average particles. This statement and the physical rationale put forward in this subsection are strongly supported by the particle-resolving numerical simulations and thought experiment presented in the next section.

3.1. Particle-resolving simulations

We use the same numerical model and code as Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017), which have undergone thorough testing and calibration in previous studies, including grid convergence tests and others (Biegert Reference Biegert2018; Yao et al. Reference Yao, Biegert, Vowinckel, Köllner, Meiburg, Balachandar, Criddle and Fringer2022; Zhu et al. Reference Zhu, He, Zhao, Vowinckel and Meiburg2022). The model is based on the immersed boundary method after Uhlmann (Reference Uhlmann2005), where each particle $p$ ’s interior is replaced by pseudofluid, while a Lagrangian surface force density $\boldsymbol{\tilde f}^p_{\mathit{IBM}}$ enforces no-slip boundary conditions at its surface $\mathbb{S}^p$ . This surface is modelled as a very thin volume layer, the so-called Lagrangian layer (Tschisgale et al. Reference Tschisgale, Kempe and Fröhlich2017), extending minimally both into the interior and exterior of the particle $p$ . A smooth volumetric force density derived from $\boldsymbol{\tilde f}^p_{\mathit{IBM}}$ also enters the Navier–Stokes equations as an additional body force term. These governing equations are then numerically integrated, on a uniform rectangular grid of length $\varDelta =0.1R$ , using a third-order low-storage Runge–Kutta scheme for time advancement and a finite-difference approach for spatial discretization. The pressure field is solved using a direct solver based on the fast Fourier transform. To integrate the equations governing the particles’ motion, the Runge–Kutta scheme subdivides the three-step fluid integration procedure into a total of $15$ substeps per fluid time step. The hydrodynamic force acting on a particle $p$ is calculated as

(3.1) \begin{align} \boldsymbol{F}^p_{\mathit{h}}={\mathrm{d}}_t\int _{\mathbb{V}^p}\rho _{\mathit{f}}\boldsymbol{u}{\mathrm{d}} V-\int _{\mathbb{S}^p}\boldsymbol{\tilde f}^p_{\mathit{IBM}}{\mathrm{d}} S-\rho _{\mathit{f}}\boldsymbol{g}V^p, \end{align}

which exploits that the instantaneous, local fluid velocity $\boldsymbol{u}$ is solved in both the exterior and interior of the particles. Furthermore, contact forces are calculated using the soft-sphere model proposed by Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017). When the distance between the surfaces of two approaching particles becomes sufficiently small, the intervening fluid is squeezed out of the gap. Once the gap width falls below approximately $2\Delta$ , this process can no longer be resolved by the fluid grid. Therefore, a lubrication model, (9) of Biegert et al. (Reference Biegert, Vowinckel and Meiburg2017), is employed to account for the unresolved hydrodynamic interactions. This model is also applied to particles rebounding after collision, when fluid is drawn back into the gap. The lubrication force is dissipative, since it always acts in the direction opposite to the relative velocity. Besides, the adaptive collision time model proposed by Kempe & Fröhlich (Reference Kempe and Fröhlich2012) is employed to account for normal contact interactions. It stretches the collision process in time to match the time step of the fluid solver.

3.1.1. Sediment transport driven by viscous flow

We use our previously reported simulation (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026) of statistically steady, uniform sediment transport driven by viscous flow (figure 1). This system is characterized by a negligible generalized-drag force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ in the bed-normal direction (explained shortly), a property that allows direct access to the generalized-buoyancy force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ , which is otherwise difficult to isolate.

Figure 1.

Snapshot of a particle-resolving simulation of statistically steady, uniform sediment transport driven by viscous flow.

The gravitational acceleration $\boldsymbol{g}$ is driving both the liquid flow and granular particles down a small slope of $3^\circ$ . Physical quantities are measured in natural units of $\rho _{\mathit{f}}$ , $(\rho _{\mathit{s}}/\rho _{\mathit{f}}-1)|\boldsymbol{g}|$ and $2R$ . The simulation was conducted for $\rho _{\mathit{s}}=2.65$ and a high viscosity $\eta _{\mathit{f}}=0.36$ , yielding a weakly inertial flow regime with a bulk flow Reynolds number $\textit{Re} \equiv \rho _{\mathit{f}}\max (u_{\mathit{f}x})H/\eta _{\mathit{f}}\approx 22.4$ , where $H$ is the extension of the vertical domain. The flow direction is $x$ , the direction normal to the bed is $z$ (henceforth called the ‘vertical’ direction for simplicity), and the lateral direction is $y$ . The system dimensions are $(L_x,L_y,H)=(5,3,8)$ , and $50$ identical, spherical particles are considered, forming a bed of approximately four particle layers, with the top layer being the most mobile (figure 1). Here $L_x$ and $L_y$ denote the length and width, respectively, of the $(x,y)$ -domain periodic in both $x$ and $y$ , while no-slip and free-slip boundary conditions are imposed at the bottom ( $z=0$ ) and top ( $z=H$ ) walls, respectively, of the domain to mimic a rigid and free surface.

The fluid’s and particles’ equation of motion are averaged using the following averaging procedure:

(3.2) \begin{align} \langle A\rangle [z_n]\equiv \frac {1}{N_TL_xL_y\Delta z}\sum _{k=1}^{N_T}\int _{z_n-\Delta z/2}^{z_n+\Delta z/2}\int _0^{L_y}\int _0^{L_x}A[\boldsymbol{x},t_k]\,{\mathrm{d}} x{\mathrm{d}} y{\mathrm{d}} z, \end{align}

which satisfies the rules in (2.4). Due to statistical bed-tangential homogeneity ( $\partial _x=\partial _y=0$ ), the averaging takes place over the entire periodic bed-tangential domain, $x\in [0,L_x)$ and $y\in [0,L_y)$ , and, due to statistical steadiness, over a sufficient number $N_T=268$ of timewise well-separated instants $t_k$ to ensure $\partial _t\approx 0$ . Furthermore, $\Delta z=0.1$ is the constant extent of the intervals with centres $z_n$ that discretize the bed-normal domain $z\in (0,H)$ . It is twice the simulation’s grid length $\varDelta =0.05$ and sufficiently small to ensure that (3.2) approximately satisfies the Reynolds rule (2.5). By contrast, since the macroscopic flow can change very rapidly over the particle diameter $2R$ in the $z$ -direction, a value of $\Delta z$ or the coarse-graining width of the order of $2R$ , as used in some previous studies (e.g. Fry et al. Reference Fry, Lacaze, Bonometti, Elyakime and Charru2024), would be too large. For reasons described by Pähtz et al. (Reference Pähtz, Chen, Zhu, Tholen and He2026), it is also undesirable to use coarse-graining kernels without a compact support, such as the often-used Gaussian kernels (Goldhirsch Reference Goldhirsch2010).

Due to $\partial _x=\partial _y=\partial _t=0$ , the average fluid phase momentum balance (2.8) in the $x$ - and $z$ -directions, respectively, becomes

(3.3) \begin{align} \frac {{\mathrm{d}}}{{\mathrm{d}} z}\big(\sigma ^{\,\mathit{f}}_{zx}+\sigma ^{\,\mathit{f}}_{\textit{Re}zx}\big)&=-\beta _{\mathit{f}}\rho _{\mathit{f}}g_x+\beta _{\mathit{s}}\langle f_{\mathit{h}x}\rangle ^{\mathit{s}},\\[-12pt]\nonumber \end{align}

(3.4) \begin{align} \frac {{\mathrm{d}}}{{\mathrm{d}} z}\big(\sigma ^{\,\mathit{f}}_{zz}+\sigma ^{\,\mathit{f}}_{\textit{Re}zz}\big)&=-\beta _{\mathit{f}}\rho _{\mathit{f}}g_z+\beta _{\mathit{s}}\langle f_{\mathit{h}z}\rangle ^{\mathit{s}}. \end{align}

Almost all coarse-grained quantities therein are computed directly as described by Pähtz et al. (Reference Pähtz, Chen, Zhu, Tholen and He2026). The only exception is $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ within $\boldsymbol{\sigma ^{\,\mathit{f}}}\equiv \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ , which is calculated indirectly from the other quantities via integrating the right-hand sides of (3.3) and (3.4) to minimize potential numerical inaccuracies, since its direct calculation is the least accurately determinable term for the immersed boundary method (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026).

Figure 2 shows the vertical profiles of the dispersed-phase volume fraction $\beta _{\mathit{s}}$ , fluid-phase-averaged streamwise flow velocity $u_{\mathit{f}x}$ , and dispersed-phase-averaged streamwise flow velocity $u_{\mathit{s}x}$ .

Figure 2.

Vertical profiles of (a) the dispersed-phase volume fraction $\beta _{\mathit{s}}$ and (b) the fluid-phase-averaged streamwise flow velocity $u_{\mathit{f}x}$ and dispersed-phase-averaged streamwise flow velocity $u_{\mathit{s}x}$ , obtained from a particle-resolving simulation of statistically steady, uniform viscous sediment transport.

All three kinds of profiles are qualitatively consistent with expectations from existing sediment transport simulations (e.g. Maurin et al. Reference Maurin, Chauchat, Chareyre and Frey2015; Pähtz & Durán Reference Pähtz and Durán2018) and highly resolved experiments (Ni & Capart Reference Ni and Capart2018). In particular, the profiles of $\beta _{\mathit{s}}$ indicate layering within the sediment bed and a bed-transport-layer interface at $z\approx 3$ (figure 2 a).

Figures 3 and 4 show the vertical profiles of the total fluid phase shear ( $\sigma ^{\,\mathit{f}}_{zx}+\sigma ^{\,\mathit{f}}_{\textit{Re}zx}$ ) and normal ( $\sigma ^{\,\mathit{f}}_{zz}+\sigma ^{\,\mathit{f}}_{\textit{Re}zz}$ ) stresses, respectively, and their constituting contributions based on $\boldsymbol{\sigma ^{\,\mathit{f}}}\equiv \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ .

Figure 3.

(a) Fluid phase shear stress contributions obtained from a particle-resolving simulation of statistically steady, uniform viscous sediment transport. (b) The sum of all contributions results in the total fluid phase shear stress $\sigma ^{\,\mathit{f}}_{zx}+\sigma _{\textit{Re}zx}^{\,\mathit{f}}=\beta _{\mathit{f}}\langle \sigma _{zx}\rangle ^{\,\mathit{f}}+\sigma _{\mathit{s}zx}^{\,\mathit{f}}+\sigma _{\textit{Re}zx}^{\,\mathit{f}}$ .

Figure 4.

(a) Fluid phase normal stress contributions obtained from a particle-resolving simulation of statistically steady, uniform viscous sediment transport. (b) The sum of all contributions results in the total fluid phase normal stress $\sigma ^{\,\mathit{f}}_{zz}+\sigma _{\textit{Re}zz}^{\,\mathit{f}}=\beta _{\mathit{f}}\langle \sigma _{zz}\rangle ^{\,\mathit{f}}+\sigma _{\mathit{s}zz}^{\,\mathit{f}}+\sigma _{\textit{Re}zz}^{\,\mathit{f}}$ , which exhibits a nearly hydrostatic profile, $\sigma ^{\,\mathit{f}}_{zz}+\sigma _{\textit{Re}zz}^{\,\mathit{f}}\approx \rho _{\mathit{f}}(H-z)g_z$ .

In particular, it can be seen that the Reynolds stress $\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}$ is negligible (figures 3 a and 4 a), as expected, and that the vertical profiles of $\sigma ^{\,\mathit{f}}_{zz}+\sigma ^{\,\mathit{f}}_{\textit{Re}zz}$ are approximately hydrostatic (figure 4 b), $\sigma ^{\,\mathit{f}}_{zz}+\sigma ^{\,\mathit{f}}_{\textit{Re}zz}\approx \rho _{\mathit{f}}(H-z)g_z$ , as was also reported by Fry et al. (Reference Fry, Lacaze, Bonometti, Elyakime and Charru2024) for their simulations of viscous sediment transport.

Due to the weakly inertial flow conditions (low Reynolds and Stokes numbers), the instantaneous hydrodynamic forces acting on the transported particles are approximately given by the background-flow (buoyancy) and drag forces in the MRG equation (1.1), and lubrication forces due to fluid-mediated particle interactions. However, as already mentioned, short-range lubrication forces are not resolved by the numerical model, but explicitly accounted for via contact force corrections, i.e. they contribute to the contact stress $\boldsymbol{\sigma ^{\mathit{s}}}$ rather than the hydrodynamic force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{h}}\rangle ^{\mathit{s}}$ in (2.8) and (2.9). Therefore, in the vertical direction, where $u_{\mathit{f}z}=u_{\mathit{s}z}=0$ due to the steady state condition, drag forces on transported particles cancel out on average and only the buoyancy force $F_{\mathit{B}z}$ remains, which is approximately equal to $-(4/3)\pi R^3\rho _{\mathit{f}}g_z$ (Maxey & Riley Reference Maxey and Riley1983). Hence, we expect

(3.5) \begin{align} \beta _{\mathit{s}}\langle f_{\mathit{h}z}\rangle ^{\mathit{s}}\approx \beta _{\mathit{s}}\langle f_{\mathit{B}z}\rangle ^{\mathit{s}}\approx -\beta _{\mathit{s}}\rho _{\mathit{f}}g_z. \end{align}

In figure 5, we use this expectation to test the buoyancy closures by Zhang et al. (Reference Zhang, VanderHeyden, Zou and Padial-Collins2007), (2.16), Ishii (Reference Ishii1975), (2.18) and Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), against the simulation data.

Figure 5.

(a) Test of the buoyancy closures by Zhang et al. (Reference Zhang, VanderHeyden, Zou and Padial-Collins2007), (2.16), Ishii (Reference Ishii1975), (2.18), and Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), against data from a particle-resolving simulation of statistically steady, uniform viscous sediment transport. If lift forces are not significant, as indicated by $\beta _{\mathit{s}}\langle f_{\mathit{h}z}\rangle ^{\mathit{s}}/(-\beta _{\mathit{s}}\rho _{\mathit{f}}g_z)\approx 1$ in (b), one expects $\beta _{\mathit{s}}\langle f_{\mathit{B}z}\rangle ^{\mathit{s}}/(\beta _{\mathit{s}}\langle f_{\mathit{h}z}\rangle ^{\mathit{s}})\approx 1$ in (a). This is, indeed, the case for the closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), but not for the other closures. The peaks around $z=4$ likely have a numerical origin due to insufficient data for statistical averaging when $\beta _{\mathit{s}}$ is very small (figure 2 a).

It can be seen that these data clearly support the latter closure. It means that the full effective mechanical fluid phase stress $\boldsymbol{\sigma ^{\,\mathit{f}}}$ , rather than only the fluid-phase-averaged stress $\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ , is the background flow stress responsible for buoyancy. This result is further confirmed by the thought experiment in § 3.2. Note that these simulation data are unable to distinguish between the closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), and the Jackson closure, (2.20), due to $\sigma ^{\,\mathit{f}}_{\textit{Re}zz}\approx 0$ .

3.1.2. ‘Horizontal settling’

This section presents original results from new simulations not reported before. We consider a single spherical particle located at the centre of a simulation box with dimensions $(L_x,L_y,H)=(20R,20R,20R)$ . The box’s boundary conditions are as in § 3.1.1: periodic in $x$ and $y$ , no-slip at $z=0$ and free-slip at $z=H$ . However, there is no longer gravity ( $\boldsymbol{g}=0$ ), but the incompressible, Newtonian fluid flow is driven by a constant negative streamwise pressure gradient, $\rho _{\mathit{f}}a\equiv -\partial _xP\gt 0$ . Physical quantities are now measured in natural units of $\rho _{\mathit{f}}$ , $a$ and $2R$ , that is, the constant acceleration $a$ replaces the acceleration scale $(\rho _{\mathit{s}}/\rho _{\mathit{f}}-1)|\boldsymbol{g}|$ used in § 3.1.1. In addition to hydrodynamic forces, the particle is subjected to a vertical force that always exactly cancels the vertical hydrodynamic force, in order to keep the vertical location constant at $z=5$ and its vertical centre-of-mass velocity constant at $v_{{\uparrow }z}=0$ . Otherwise, for turbulent conditions, the particle would migrate towards the bottom wall, making it much more challenging to obtain unambiguous results. Since the vertical force does not do work, it represents a purely kinematic constraint that does not affect the flow forces acting on the particle, analogous to kinematic constraints in analytical mechanics. It could, for example, be established through giving the particle a charge and imposing an electric field that instantly recognizes and adapts to the hydrodynamics. Note that instant signalling and action are consistent with Newton’s physics, which fluid dynamics is based on.

The above system has been designed to mimic the vertical settling of a particle in a quiescent fluid, but in the horizontal direction (figure 6).

Figure 6.

Sketch of the horizontal settling scenario.

While for the vertical settling scenario, the generalized-buoyancy force emerges due to vertical normal stress gradients, for this ‘horizontal settling’ scenario, it emerges due to vertical shear stress gradients and the flow-driving horizontal pressure gradient. However, as we argued in § 2.4 and as predicted by the closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), Reynolds shear stress gradients, ${\mathrm{d}}_z\sigma ^{\,\mathit{f}}_{\textit{Re}zx}$ , should not contribute to buoyancy. Therefore, for turbulent flow, where the Reynolds shear stress $\sigma ^{\,\mathit{f}}_{\textit{Re}zx}$ dominates the effective mechanical shear stress in the average fluid phase momentum balance,

(3.6) \begin{align} \rho _{\mathit{f}}a=-\frac {{\mathrm{d}}}{{\mathrm{d}} z}\big(\sigma ^{\,\mathit{f}}_{zx}+\sigma ^{\,\mathit{f}}_{\textit{Re}zx}\big), \end{align}

we expect that the shear stress gradient force is negligible. Then, in the steady state, the only two significant forces acting on the particle on average are the force from the constant pressure gradient, $(4/3)\pi R^3\rho _{\mathit{f}}a$ , and the drag force. Hence, the particle should lead the flow by approximately the settling velocity $v_{\mathit{settle}}$ , which can be calculated by (9) of Dietrich (Reference Dietrich1982),

(3.7) \begin{align} \log _{10}v_{\mathit{settle}} & =-3.76715+1.92944(\log _{10}\mathit{Ar})-0.09815(\log _{10}\mathit{Ar})^2 \nonumber\\[3pt] & \quad -0.00575(\log _{10}\mathit{Ar})^3+0.00056(\log _{10}\mathit{Ar})^4, \end{align}

where $\mathit{Ar}\equiv \rho _{\mathit{f}}^2a(2R)^3/\eta _{\mathit{f}}^2$ (i.e. $\mathit{Ar}=1/\eta _{\mathit{f}}^2$ in natural units) is the Archimedes number. For the averaged velocity fields evaluated at the particle centre location, $z=5$ , this means

(3.8) \begin{align} \text{(turbulent flow)}\qquad (u_{\mathit{f}x}-u_{\mathit{s}x})|_{z=5}\approx -v_{\mathit{settle}}. \end{align}

By contrast, if the Reynolds stress actually contributed to buoyancy, as predicted by the Jackson closure, (2.20), the shear stress gradient force would compensate the pressure gradient force $(4/3)\pi R^3\rho _{\mathit{f}}a$ due to (3.6). Then, only the drag force would act on the particle, and we would therefore expect $u_{\mathit{f}x}-u_{\mathit{s}x}\approx 0$ at $z=5$ . Furthermore, regardless of the chosen buoyancy closure, for fully laminar flow, the shear stress gradient force always exactly compensates the pressure gradient force. For such conditions, not only the bulk flow Reynolds number $\textit{Re}$ but also the particle Reynolds numbers $\textit{Re} _U$ and $\textit{Re} _W$ are small, meaning that the assumptions behind the MRG equation (1.1) are approximately satisfied. Then, the counteraction of the shear stress gradient force $F_{\tau x}$ is further enhanced by the Faxén correction to the drag force by a value that is equal to $(3/4)F_{\tau x}$ . This is because the surface integral $\int _{\mathbb{S}}(\boldsymbol{\tilde u}-\boldsymbol{v}_{\uparrow} ){\mathrm{d}} S$ in (1.1) simplifies to $4\pi R^2(\boldsymbol{\tilde u}[\boldsymbol{x}_{\mathit{o}}]-\boldsymbol{v}_{\uparrow} +R^2\Delta \boldsymbol{\tilde u}[\boldsymbol{x}_{\mathit{o}}]/6)$ , since the steady Stokes equations imply $\Delta ^k\boldsymbol{\tilde u}=0$ for any integer $k\gt 1$ (Kim & Karrila Reference Kim and Karrila1991). For fully laminar flow, we therefore expect

(3.9) \begin{align} \text{(laminar flow)}\qquad (u_{\mathit{f}x}-u_{\mathit{s}x})|_{z=5}=+(3/4)v_{\mathit{settle}}. \end{align}

The three strikingly different expectations, depending on the flow conditions and chosen buoyancy closure, for the outcome of $(u_{\mathit{f}x}-u_{\mathit{s}x})|_{z=5}$ constitute an excellent environment for testing the different closures with our particle-resolving simulations. We simulate both a viscous scenario, $\mathit{Ar}=1$ , and a turbulent scenario, $\mathit{Ar}=10^4$ . The corresponding bulk flow Reynolds numbers are $\textit{Re} \approx 499$ and $\textit{Re} \approx 78\,092$ , respectively. To allow a quick adaptation of the particle to the flow velocity, we choose a low particle density, $\rho _{\mathit{s}}=0.5$ , which is close to the smallest value that does not lead to numerical issues. As in § 3.1.1, the simulations’ grid length $\varDelta =0.05$ . Once the simulations have reached a steady state, the data are averaged as in (3.2) over the entire periodic bed-tangential domain, $x\in [0,L_x)$ and $y\in [0,L_y)$ , and over a sufficient number $N_T$ of timewise well-separated instants $t_k$ , $N_T=214$ for the viscous and $N_T=2831$ for the turbulent case.

Figure 7 confirms that the effective mechanical fluid phase shear stress $\sigma ^{\,\mathit{f}}_{zx}$ dominates for the viscous case, whereas the Reynolds shear stress $\sigma ^{\,\mathit{f}}_{\textit{Re}zx}$ dominates for the turbulent case. Furthermore, figure 8 shows the time series of the spatially averaged fluid phase ( $u_{\mathit{f}x,t}$ ) and dispersed-phase ( $u_{\mathit{s}x,t}$ ) streamwise flow velocities, and their difference $u_{\mathit{f}x,t}-u_{\mathit{s}x,t}$ , all evaluated at the particle centre elevation $z=5$ . It can be seen that the mean velocity difference $(u_{\mathit{f}x}-u_{\mathit{s}x})|_{z=5}$ (cyan lines in figure 8 c, f) is positive ( $\approx 1.1v_{\mathit{settle}}$ ) in the viscous case and negative ( $\approx -1.2v_{\mathit{settle}}$ ) in the turbulent case. Since these values are not far off the predictions $(3/4)v_{\mathit{settle}}$ and $-v_{\mathit{settle}}$ by (3.9) and (3.8), respectively, they strongly support the buoyancy closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), over the Jackson closure, (2.20). In addition, the vertical profiles of $u_{\mathit{f}x}$ and $u_{\mathit{s}x}$ reveal that $u_{\mathit{f}x}$ remains smaller than $u_{\mathit{s}x}$ for all locations $z\in (4.5,5.5)$ covered by the particle in the turbulent case (figure 9 b). This would not be physically possible if a Reynolds shear stress gradient force compensated the pressure gradient force. Finally, figure 10 confirms that the findings for the turbulent case are neither the result of correlations between the turbulent fluid velocity fluctuations and drag forces – this possibility is ruled out by the fact that the distribution of $u_{\mathit{f}x,t}|_{z=5}$ is nearly symmetric relative to its mean $u_{\mathit{f}x}|_{z=5}$ (figure 10 a) – nor caused by insufficient averaging (figure 10 b). Put together, figures 8, 9 and 10 falsify the Jackson closure, (2.20).

Figure 7.

Fluid phase shear stress contributions obtained from particle-resolving simulations of the horizontal settling scenario (figure 6) for (a) viscous flow ( $\mathit{Ar}=1$ , $\textit{Re} \approx 499$ ) and (b) turbulent flow ( $\mathit{Ar}=10^4$ , $\textit{Re} \approx 78\,092$ ).

Figure 8.

Time series of the spatially averaged fluid phase ( $u_{\mathit{f}x,t}$ ) and dispersed-phase ( $u_{\mathit{s}x,t}$ ) streamwise flow velocities, and their difference $u_{\mathit{f}x,t}-u_{\mathit{s}x,t}$ , all evaluated at the particle centre elevation $z=5$ . The cyan lines indicate the respective mean values. The data are from particle-resolving simulations of the horizontal settling scenario (figure 6) for (a–c) viscous flow ( $\mathit{Ar}=1$ , $\textit{Re} \approx 499$ ) and (d–f) turbulent flow ( $\mathit{Ar}=10^4$ , $\textit{Re} \approx 78\,092$ ).

Figure 9.

Vertical profiles of the fluid-phase-averaged streamwise flow velocity $u_{\mathit{f}x}$ and the dispersed-phase-averaged streamwise flow velocity $u_{\mathit{s}x}$ , obtained from particle-resolving simulations of the horizontal settling scenario (figure 6) for (a) viscous flow ( $\mathit{Ar}=1$ , $\textit{Re} \approx 499$ ) and (b) turbulent flow ( $\mathit{Ar}=10^4$ , $\textit{Re} \approx 78\,092$ ).

3.2. Thought experiment: dilute Stokesian suspension

We consider a mixture of an incompressible Newtonian fluid and identical spherical particles at vanishingly small particle Reynolds numbers and Stokes number, $\textit{Re} _U\sim 0$ , $\textit{Re} _W\sim 0$ and $\mathit{St}\equiv \rho _{\mathit{s}}R^2/(T\eta _{\mathit{f}})\sim 0$ , where $T$ is the characteristic time over which the disturbance flow velocity changes. It means that transients in the Stokes disturbance flow for any given particle can also be neglected. Furthermore, the dispersed-phase volume fraction $\beta _{\mathit{s}}$ shall be non-zero but sufficiently low to ensure that fluid-mediated particle–particle interactions are negligible, while the sphere radius $R$ shall be sufficiently smaller than the characteristic length $L$ over which the macroscopic flow changes. Under these assumptions, the following average fluid phase momentum balance can be derived ((2.46) of Jackson Reference Jackson2000):

(3.10) \begin{align} \rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\left \{-\langle P\rangle ^{\,\mathit{f}}\boldsymbol{1}+\eta _{\mathit{f}}\big(\boldsymbol{\nabla }\langle \boldsymbol{u}\rangle +(\boldsymbol{\nabla }\langle \boldsymbol{u}\rangle )^{\textrm{T}}\big)+\frac {5\eta _{\mathit{f}}\beta _{\mathit{s}}}{2}\big(\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}}+(\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}})^{\textrm{T}}\big)\right . \nonumber\\[3pt] & \quad +\left .3\eta _{\mathit{f}}\beta _{\mathit{s}}\left (\frac {1}{2}\big(\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}}-\left (\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}}\right )^{\textrm{T}}\big)-\boldsymbol{\epsilon }\boldsymbol{\cdot }\langle \boldsymbol{\omega }\rangle ^{\mathit{s}}\right )-\boldsymbol{\nabla }\left [\frac {3\eta _{\mathit{f}}\beta _{\mathit{s}}}{4}\left (\boldsymbol{u}_{\mathit{f}}-\boldsymbol{u}_{\mathit{s}}\right )\right ]\right \} \nonumber\\[3pt] & \quad -\frac {9\eta _{\mathit{f}}\beta _{\mathit{s}}}{2R^2}\left (\boldsymbol{u}_{\mathit{f}}-\boldsymbol{u}_{\mathit{s}}+\frac {R^2}{6}\Delta \boldsymbol{u}_{\mathit{f}}\right )+\rho _{\mathit{f}}\boldsymbol{g}, \end{align}

where $\varDelta =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\nabla }$ is the Laplace operator, $\boldsymbol{\epsilon }$ the total antisymmetric tensor, $P$ the instantaneous, local fluid pressure and $\boldsymbol{\omega }\equiv \sum _pX^p\boldsymbol{\omega }^p$ the particle-angular velocity field. Note that, for the considered scenario, $\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}=0$ due to $\textit{Re} _W\sim 0$ and $\mathit{St}\sim 0$ (Jackson Reference Jackson1997) and $\langle \boldsymbol{u}\rangle =\beta _{\mathit{f}}\boldsymbol{u}_{\mathit{f}}+\beta _{\mathit{s}}\boldsymbol{u}_{\mathit{s}}$ due to constant microscopic fluid and particle densities.

Figure 10.

(a) Histogram of the spatially averaged fluid phase velocity $u_{\mathit{f}x,t}$ relative to its mean $u_{\mathit{f}x}$ , evaluated at the particle centre elevation $z=5$ . (b) Temporal evolution of the velocity difference averaged up to time $t$ , $(u_{\mathit{f}x}^t-u_{\mathit{s}x}^t)|_{z=5}$ and its standard error. The data are from particle-resolving simulations of the horizontal settling scenario (figure 6) for turbulent flow ( $\mathit{Ar}=10^4$ , $\textit{Re} \approx 78\,092$ ).

The question is now how to interpret the different terms in (3.10). Due to $\textit{Re} _U\sim 0$ , $\textit{Re} _W\sim 0$ and $\mathit{St}\sim 0$ , the disturbance flow is described by the steady Stokes equations. Hence, the instantaneous disturbance or generalized-drag force acting on any given particle centred at $\boldsymbol{x}_{\mathit{o}}$ is described by only the drag force and its Faxén correction in the MRG equation (1.1). In particular, the surface integral $\int _{\mathbb{S}}(\boldsymbol{\tilde u}-\boldsymbol{v}_{\uparrow} ){\mathrm{d}} S$ in (1.1) simplifies to $4\pi R^2(\boldsymbol{\tilde u}[\boldsymbol{x}_{\mathit{o}}]-\boldsymbol{v}_{\uparrow} +R^2\varDelta \boldsymbol{\tilde u}[\boldsymbol{x}_{\mathit{o}}]/6)$ , since the steady Stokes equations imply $\Delta ^k\boldsymbol{\tilde u}=0$ for any integer $k\gt 1$ (Kim & Karrila Reference Kim and Karrila1991). Therefore, identifying $\boldsymbol{u}_{\mathit{f}}$ and $\boldsymbol{u}_{\mathit{s}}$ as the analogues to $\boldsymbol{\tilde u}[\boldsymbol{x}_{\mathit{o}}]$ and $\boldsymbol{v}_{\uparrow}$ in (1.1), respectively, the first term on the third line of (3.10) is recognized as the negative sum of the average drag force and its Faxén correction per unit volume of the pure fluid. However, the actual generalized-drag force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ is the sum of the average drag force and its Faxén correction per unit volume of the fluid–particle mixture, not of the pure fluid. This is because fluid phase momentum balances in two-fluid models refer to the mixture volume. To translate pure-fluid-based volumetric densities into mixture-based volumetric densities, one has to multiply the former with the fluid volume fraction $\beta _{\mathit{f}}$ . Hence, we conclude

(3.11) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}=\frac {9\eta _{\mathit{f}}\beta _{\mathit{s}}\beta _{\mathit{f}}}{2R^2}\left (\boldsymbol{u}_{\mathit{f}}-\boldsymbol{u}_{\mathit{s}}+\frac {R^2}{6}\Delta \boldsymbol{u}_{\mathit{f}}\right )\!. \end{align}

The remaining part $(\beta _{\mathit{s}}^2/\beta _{\mathit{f}})\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ is therefore attributed to the generalized-buoyancy force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ alongside the usual gravitational and inertial terms:

(3.12) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}=\beta _{\mathit{s}}\left ({\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}\right )\!. \end{align}

This somewhat odd appearance of drag terms in $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ was exactly what led to the debate between Clift et al. (Reference Clift, Seville, Moore and Chavarie1987) and others that we briefly discussed in § 1.

The remaining terms in (3.10) are the gravitational body force density $\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}$ and a divergence of two distinct stress tensors: the fluid-phase-averaged stress tensor (Jackson Reference Jackson1997),

(3.13) \begin{align} \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}=-\beta _{\mathit{f}}\langle P\rangle ^{\,\mathit{f}}\boldsymbol{1}+\eta _{\mathit{f}}\Big(\boldsymbol{\nabla }\langle \boldsymbol{u}\rangle +(\boldsymbol{\nabla }\langle \boldsymbol{u}\rangle )^{\textrm{T}}\Big), \end{align}

and a pseudostress tensor $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ that emerges due to fluid–particle interactions (Jackson Reference Jackson1997),

(3.14) \begin{align} \boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}=&-\beta _{\mathit{s}}\langle P\rangle ^{\,\mathit{f}}\boldsymbol{1}+\frac {5\eta _{\mathit{f}}\beta _{\mathit{s}}}{2}\Big(\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}}+\big(\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}}\big)^{\textrm{T}}\Big) \nonumber\\ &+3\eta _{\mathit{f}}\beta _{\mathit{s}}\left (\frac {1}{2}\big(\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}}-\big(\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}}\big)^{\textrm{T}}\big)-\boldsymbol{\epsilon }\boldsymbol{\cdot }\langle \boldsymbol{\omega }\rangle ^{\mathit{s}}\right )-\boldsymbol{\nabla }\left [\frac {3\eta _{\mathit{f}}\beta _{\mathit{s}}}{4}\left (\boldsymbol{u}_{\mathit{f}}-\boldsymbol{u}_{\mathit{s}}\right )\right ]\!. \end{align}

Using (3.11)–(3.14) and $\boldsymbol{\sigma ^{\,\mathit{f}}}\equiv \beta _{\mathit{f}}\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s}}}$ , we can rewrite (3.10) in either of the following two forms:

(3.15) \begin{align} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}, \\[-12pt]\nonumber \end{align}

(3.16) \begin{align} \rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}+\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}. \end{align}

Multiplying (3.16) with $\beta _{\mathit{f}}$ and subtracting the result from (3.15) then yields the buoyancy closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), here equivalent to the Jackson closure, (2.20), as $\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}=0$ for the considered scenario.

3.3. Conclusion from the particle-resolving simulations and thought experiment

The results presented in §§ 3.1 and 3.2 strongly support the buoyancy closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), over the closure by Zhang et al. (Reference Zhang, VanderHeyden, Zou and Padial-Collins2007), (2.16) and the Jackson closure, (2.20). However, as we argued in § 2.4, (2.19) constitutes only the small-particle approximation of what we consider to be the physically appropriate closure. As such, it leads to an inconsistency for some fluid–particle systems at complete rest (§ 2.2). It is therefore generalized to arbitrarily large particles in the next section.

4.1. Operator notation

The operators discussed in this section map an arbitrary macroscopic tensor fields $\boldsymbol{A}$ (i.e. fields that satisfy $\langle \boldsymbol{A}\rangle =\boldsymbol{A}$ ) in space and time, $\mathbb{R}^3\times \mathbb{R}$ , to new macroscopic tensor fields in $\mathbb{R}^3\times \mathbb{R}$ of the same rank. Their actions are the same for each of the $1$ (scalar), $3$ (vector) or $9$ (rank- $2$ tensor) components of $\boldsymbol{A}$ .

We define the operator $\mathcal{L}$ through the action

(4.1) \begin{align} \mathcal{L}\beta _{\mathit{s}}\boldsymbol{A}\equiv \left \langle \sum \nolimits _pX^p\underline {\boldsymbol{A}}_{\mathbb{V}^p}\right \rangle =\left \langle X_{\mathit{s}}\sum \nolimits _pX^p\underline {\boldsymbol{A}}_{\mathbb{V}^p}\right \rangle =\beta _{\mathit{s}}\left \langle \sum \nolimits _pX^p\underline {\boldsymbol{A}}_{\mathbb{V}^p}\right \rangle ^{\mathit{s}}, \end{align}

where we introduced the notation $\underline {\boldsymbol{\cdot }}_{\mathbb{V}^p}\equiv ({1}/{V^p})\int _{\mathbb{V}^p}\boldsymbol{\cdot }{\mathrm{d}} V$ and used $X_{\mathit{s}}X^p=X^p$ and (2.6). Equation (4.1) means that $\mathcal{L}$ first averages $\boldsymbol{A}$ over a given particle $p$ ’s domain $\mathbb{V}^p$ . Then, it assigns the resulting average to each point of the same domain, $\boldsymbol{x}\in \mathbb{V}^p$ , but to nowhere else. It does so for all particles $p$ , superposes all these new fields and macroscopically averages this superposition. The space–time dependence of the resulting field emerges from the generalized functions $X^p[\boldsymbol{x},t]$ within the averaging procedure $\langle \boldsymbol{\cdot }\rangle$ .

The operator $\mathcal{L}$ and the identity operator $\mathcal{I}$ appear in the identity (Supplementary material of Pähtz et al. (Reference Pähtz, Chen, Zhu, Tholen and He2026))

(4.2) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }[\boldsymbol{A}]_{\mathit{s}}^\prime =(\mathcal{L}-\mathcal{I})\beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{A}, \end{align}

where $[\boldsymbol{\cdot }]_{\mathit{s}}^\prime$ is another operator. Its detailed micromechanical action is given by (A7), but it is not relevant for what comes later. It is only relevant that $[\boldsymbol{\cdot }]_{\mathit{s}}^\prime$ exists. Note that, in the small-particle approximation, $[\boldsymbol{\cdot }]_{\mathit{s}}^\prime \approx 0$ and $\mathcal{L}\approx \mathcal{I}$ .

For the coming analysis, we introduce the convention that operators written using capital letters, such as $\mathcal{L}$ and $\mathcal{I}$ , shall act on all terms right of them. The exception is when both capital-letter operator and operated field are surrounded by parentheses; for example, $(\mathcal{L}\beta _{\mathit{s}})\boldsymbol{A}=\beta _{\mathit{s}}\boldsymbol{A}$ due to $\mathcal{L}\beta _{\mathit{s}}=\beta _{\mathit{s}}$ , whereas $(\mathcal{L}-\mathcal{I})\beta _{\mathit{s}}\boldsymbol{A}$ should be read as $\mathcal{L}\beta _{\mathit{s}}\boldsymbol{A}-\beta _{\mathit{s}}\boldsymbol{A}$ . Furthermore, we introduce the placeholder symbol ‘ $\underline {\,\,}$ ’, which shall indicate that the operator action is not done even when surrounded by parentheses; for example, $(\mathcal{L}\beta _{\mathit{s}}\underline {\,\,})\boldsymbol{A}=\mathcal{L}\beta _{\mathit{s}}\boldsymbol{A}$ . Objects such as $\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}$ can therefore be interpreted as operators themselves. Moreover, when the operation $\boldsymbol{\cdot }^{-1}$ is applied on an operator $\mathcal{A}$ , it shall refer to its inverse, defined through the property $\mathcal{A}^{-1}\mathcal{A}=\mathcal{A}\mathcal{A}^{-1}=\mathcal{I}$ , while $\mathcal{A}^k$ and $[\boldsymbol{\cdot }]_{\mathit{s}}^{\prime k}$ shall mean $k$ applications of $\mathcal{A}$ and $[\boldsymbol{\cdot }]_{\mathit{s}}^\prime$ , respectively.

4.2. Derivation of closure

The definition (4.1) of the operator $\mathcal{L}$ allow us to rewrite (2.22) as

(4.3) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}=\beta _{\mathit{s}}\left \langle \sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{S}^p}\boldsymbol{n}^p\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}{\mathrm{d}} S\right \rangle ^{\mathit{s}}=\left \langle \sum \nolimits _p\frac {X^p}{V^p}\int _{\mathbb{V}^p}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}{\mathrm{d}} V\right \rangle =\mathcal{L}\beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}, \end{align}

where we used $X_{\mathit{s}}X^p=X^p$ and (2.6). Then, we can restate the buoyancy criterion (2.21) in three different manners as

(4.4) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}&=\beta _{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\mathcal{L}\beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}, \\[-12pt]\nonumber \end{align}

(4.5) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\left [\boldsymbol{\hat \sigma ^{\,\mathit{f}}}\right ]_{\mathit{s}}^\prime \right )\!, \\[-12pt]\nonumber \end{align}

(4.6) \begin{align} \mathcal{L}\beta _{\mathit{s}}\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}&=\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}+\mathcal{L}\beta _{\mathit{s}}\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}, \end{align}

where (4.5) uses (4.2) and (4.6) follows from applying $\mathcal{L}\beta _{\mathit{s}}\beta _{\mathit{f}}^{-1}\underline {\,\,}$ on (4.4) using (4.3). Inverting these relations yields

(4.7) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}&=\left (\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}+\mathcal{I}\beta _{\mathit{f}}\underline {\,\,}\right )^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}=\sum _{k=0}^\infty (-1)^k\left (\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}-\mathcal{I}\beta _{\mathit{s}}\underline {\,\,}\right )^k\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}, \\[-12pt]\nonumber \end{align}

(4.8) \begin{align} \boldsymbol{\hat \sigma ^{\,\mathit{f}}}&=\left ([\boldsymbol{\cdot }]^\prime _{\mathit{s}}+\mathcal{I}\right )^{-1}\boldsymbol{\sigma ^{\,\mathit{f}}}=\sum _{k=0}^\infty (-1)^k\left [\boldsymbol{\sigma ^{\,\mathit{f}}}\right ]_{\mathit{s}}^{\prime k}, \\[-12pt]\nonumber \end{align}

(4.9) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}&=\big(\mathcal{L}\beta _{\mathit{s}}\beta _{\mathit{f}}^{-1}\underline {\,\,}+\mathcal{I}\big)^{-1}\mathcal{L}\beta _{\mathit{s}}\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}=\sum _{k=0}^\infty (-1)^k\big (\mathcal{L}\beta _{\mathit{s}}\beta _{\mathit{f}}^{-1}\underline {\,\,}\big )^{k+1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}, \end{align}

respectively, where the right-hand sides exploit the Neumann series representation of operators. Convergence of these series is often nearly guaranteed in (4.7) and (4.8) as $\beta _{\mathit{s}}$ is usually substantially smaller than unity and $[\boldsymbol{\sigma ^{\,\mathit{f}}}]_{\mathit{s}}^{\prime k}\sim \beta _{\mathit{s}}^k$ , while $\mathcal{L}-\mathcal{I}$ and $[\boldsymbol{\cdot }]^\prime _{\mathit{s}}$ typically do not increase the magnitude of a field. However, (4.9) requires at least $\beta _{\mathit{s}}\beta _{\mathit{f}}^{-1}\lt 1$ for convergence, i.e. $\beta _{\mathit{s}}\lt 0.5$ . Equation (4.8) means that a stress $\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ that satisfies the buoyancy criterion (2.21) does, indeed, exist, justifying the physical rationale laid out in § 2.4. In particular, for fluid–particle systems at complete rest, where $\boldsymbol{\sigma ^{\,\mathit{f}}}=\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}+[\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}]^\prime _{\mathit{s}}$ due to (A1) and $\boldsymbol{\sigma ^{\,\mathit{f}}_{\mathit{s^\prime }}}=0$ , it reduces to $ \boldsymbol{\hat \sigma ^{\,\mathit{f}}}=\langle \boldsymbol{\sigma }\rangle ^{\,\mathit{f}}$ , as required. In addition, (4.9) constitutes a first explicit expression (but with a not necessarily convergent series expansion) linking $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}$ to the effective mechanical stress $\boldsymbol{\sigma ^{\,\mathit{f}}}$ . Another equivalent explicit expression (with an almost always convergent series expansion) results from applying $\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}$ on (4.7) using (4.3):

(4.10) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}=\mathcal{L}\beta _{\mathit{s}}\left (\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}+\mathcal{I}\beta _{\mathit{f}}\underline {\,\,}\right )^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}. \end{align}

Yet another equivalent explicit expression is obtained from rearranging (2.23) for $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ and substituting the result in (4.3):

(4.11) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}=\mathcal{L}\beta _{\mathit{s}}\big(\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}\big ). \end{align}

This form is particularly useful for practical applications, since it involves the operator $\mathcal{L}$ only once and not in its inverted form. Yet another equivalent buoyancy expression follows from combining (2.21) and (4.4):

(4.12) \begin{align} \beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{B}}\rangle ^{\mathit{s}}=\beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\left [\boldsymbol{\hat \sigma ^{\,\mathit{f}}}\right ]_{\mathit{s}}^\prime . \end{align}

4.3. Resulting two-fluid momentum balances

Substituting the three equivalent expressions (4.10)–(4.12) in (2.8) and (2.9) yields three equivalent versions of the average two-fluid momentum balances. Equations (4.7) and (4.10) result in

(4.13) \begin{align} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\beta _{\mathit{f}}\left (\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}+\mathcal{I}\beta _{\mathit{f}}\underline {\,\,}\right )^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}, \\[-12pt]\nonumber \end{align}

(4.14) \begin{align} \beta _{\mathit{s}}\rho _{\mathit{s}}{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}&=\mathcal{L}\beta _{\mathit{s}}\left (\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}+\mathcal{I}\beta _{\mathit{f}}\underline {\,\,}\right )^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}+\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{g}+\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}. \end{align}

Equation (4.11) leads to

(4.15) \begin{align} \begin{split} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}} \\ &\quad-\mathcal{L}\beta _{\mathit{s}}\big (\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}\big ), \end{split} \\[-12pt]\nonumber \end{align}

(4.16) \begin{align} \begin{split} \beta _{\mathit{s}}\rho _{\mathit{s}}{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}+\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{g}+\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}} \\ &\quad+\mathcal{L}\beta _{\mathit{s}}\big (\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}\big ). \end{split} \end{align}

Equation (4.12) yields

(4.17) \begin{align} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\beta _{\mathit{f}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}, \\[-12pt]\nonumber \end{align}

(4.18) \begin{align} \beta _{\mathit{s}}\rho _{\mathit{s}}{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}&=\beta _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\mathit{s}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}+\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{g}+\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}, \end{align}

where (4.17) is the same as (2.23) and $\boldsymbol{\hat \sigma ^{\mathit{s}}}\equiv \boldsymbol{\sigma ^{\mathit{s}}}+[\boldsymbol{\hat \sigma ^{\,\mathit{f}}}]_{\mathit{s}}^\prime$ is a modified effective mechanical dispersed-phase stress tensor. Furthermore, (4.17) and (4.18) can be rearranged and brought into the form

(4.19) \begin{align} \begin{split} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}} \\ &\quad-\beta _{\mathit{s}}\Big (\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}\Big ), \end{split} \\[-12pt]\nonumber \end{align}

(4.20) \begin{align} \begin{split} \beta _{\mathit{s}}\rho _{\mathit{s}}{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}&=\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\mathit{s}}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}+\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{g}+\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}} \\ &\quad+\beta _{\mathit{s}}\Big (\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}\Big ). \end{split} \end{align}

The former two versions of the average two-fluid momentum balances, (4.13)–(4.16), are based on the effective mechanical stress tensors $\boldsymbol{\sigma ^{\,\mathit{f}}}$ and $\boldsymbol{\sigma ^{\mathit{s}}}$ . They are of a mathematical structure never reported before due to the explicit appearance of the operator $\mathcal{L}$ . By contrast, the latter two versions, (4.17)–(4.20), are of the same mathematical structure as those reported by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), since the effects of $\mathcal{L}$ are buried in the modified effective mechanical stress tensors $\boldsymbol{\hat \sigma ^{\,\mathit{f}}}$ and $\boldsymbol{\hat \sigma ^{\mathit{s}}}$ . Note that $\boldsymbol{\hat \sigma ^{\mathit{s}}}$ does not only differ quantitatively but also qualitatively from $\boldsymbol{\sigma ^{\mathit{s}}}$ , since, in contrast to $\boldsymbol{\sigma ^{\mathit{s}}}$ , it contains a stress contribution, $[\boldsymbol{\hat \sigma ^{\,\mathit{f}}}]_{\mathit{s}}^\prime$ , that originated in the fluid phase momentum balance.

4.4. Effect of $\mathcal{L}$ for a mixture of continuous fluid and identically sized spheres

The mathematical effect of the operator $\mathcal{L}$ can be fully understood for the simplified case of a disperse two-phase flow in which all particles are spheres of the same radius $R$ and volume $V=4\pi R^3/3$ . In this case, the action of $\mathcal{L}$ can be expressed in real and Fourier space, respectively, as (Appendix B)

(4.21) \begin{align} \mathcal{L}\beta _{\mathit{s}}\boldsymbol{A}&=\mathcal{M}\big (\mathcal{M}^{-1}\beta _{\mathit{s}}\big )\mathcal{M}\boldsymbol{A}, \\[-12pt]\nonumber \end{align}

(4.22) \begin{align} \mathcal{F}\mathcal{L}\beta _{\mathit{s}}\boldsymbol{A}&=M_K\big (\big (M_K^{-1}\mathcal{F}\beta _{\mathit{s}}\big )\ast \big (M_K\mathcal{F}\boldsymbol{A}\big )\big), \end{align}

where $\mathcal{F}$ denotes the Fourier transformation with respect to $\boldsymbol{x}$ , ‘ $\ast$ ’ the convolution, $\mathcal{M}$ is the volume averaging operator defined by (B3) or through $\mathcal{M}\boldsymbol{A}[\boldsymbol{x}]=({1}/{V})\int _{|\boldsymbol{r}|\lt R}\boldsymbol{A}[\boldsymbol{x}+\boldsymbol{r}]{\mathrm{d}} V_r$ , and $M_K$ its Fourier transform given by

(4.23) \begin{align} M_K=\sum _{j=0}^\infty \frac {3(\mathrm{i}K)^{2j}}{(2j+1)!(2j+3)}=\frac {3(\sin K-K\cos K)}{K^3}, \end{align}

with $K\equiv R|\boldsymbol{k}|$ the non-dimensionalized absolute value of the wavenumber vector $\boldsymbol{k}$ spanning the Fourier space. The effect of $\mathcal{L}$ becomes clear once (4.22) is compared with

(4.24) \begin{align} \mathcal{F}\beta _{\mathit{s}}\boldsymbol{A}=(\mathcal{F}\beta _{\mathit{s}})\ast (\mathcal{F}\boldsymbol{A}). \end{align}

It can be seen that $\mathcal{L}$ modulates the field $\beta _{\mathit{s}}\boldsymbol{A}$ in Fourier space via the function $M_K$ , which is essentially a low-pass filter that cuts off large- $K$ contributions, while approaching unity in the limit $K\to 0$ (figure 11).

Figure 11.

The Fourier transform $M_K$ of the volume averaging operator $\mathcal{M}$ essentially acts like a low-pass filter, cutting off contributions from wavenumber vectors $\boldsymbol{k}$ with $K\equiv R|\boldsymbol{k}|\gtrsim 4$ .

Note that $\mathcal{M}^{-1}\beta _{\mathit{s}}=nV$ , where $n$ is the number density, and, thus, $M_K^{-1}\mathcal{F}\beta _{\mathit{s}}=V\mathcal{F}n$ , implying $M_K^{-1}\mathcal{F}\beta _{\mathit{s}}$ in (4.22) is well-defined for values of $\boldsymbol{k}$ where $M_K=0$ due to $\mathcal{F}\beta _{\mathit{s}}[\boldsymbol{k}]=0$ .

Choosing $\boldsymbol{A}=\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ in (4.24) results in the Fourier transform of the closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), whereas the same choice in (4.22) yields the Fourier transform of its generalization (4.11). Hence, at least for identical spherical particles, the latter is essentially a low-pass-filtered version of the former closure.

5.1. Linear well-posedness of two-fluid models

Linear well-posedness is one of three criteria, next to the existence and uniqueness of a solution, required for a PDE initial-value problem to be well-posed in the sense of Hadamard (Pinchover & Rubinstein Reference Pinchover and Rubinstein2005). To analyse whether two-fluid models are linearly well-posed when employing our buoyancy closure, we consider a two-fluid model with the closure relations

(5.1) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}}&=-\boldsymbol{\nabla }P_{\mathit{f}}+\eta _{\mathit{m}}\Delta \boldsymbol{u}_{\mathit{m}},&\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\mathit{s}}}&=-P_{\mathit{s}}^\prime \boldsymbol{\nabla }\beta _{\mathit{s}},&\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}&=0,&\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}&=0,&\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}&=0, \end{align}

where $P_{\mathit{f}}$ is the effective fluid phase pressure, $P_{\mathit{s}}[\beta _{\mathit{s}}]$ the $\beta _{\mathit{s}}$ -dependent effective dispersed-phase pressure (e.g. see Johnson et al. Reference Johnson, Nott and Jackson1990; Rauter Reference Rauter2021), with the physical property $P^\prime _{\mathit{s}}\geqslant 0$ (where the prime denotes the derivative), $\eta _{\mathit{m}}[\beta _{\mathit{s}}]$ the $\beta _{\mathit{s}}$ -dependent mixture viscosity and $\boldsymbol{u}_{\mathit{m}}=\beta _{\mathit{f}}\boldsymbol{u}_{\mathit{f}}+\beta _{\mathit{s}}\boldsymbol{u}_{\mathit{s}}$ the mixture velocity. We choose the form $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau ^{\,\mathit{f}}}=\eta _{\mathit{m}}\Delta \boldsymbol{u}_{\mathit{m}}$ as viscous dissipation, since it emerges analytically for particulate Stokes flow (Jackson Reference Jackson1997) and is also used to model more complex disperse two-phase flows (e.g. Revil-Baudard & Chauchat Reference Revil-Baudard and Chauchat2013). However, the more conventional form $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau ^{\,\mathit{f}}}=\eta _{\mathit{f}}\Delta \boldsymbol{u}_{\mathit{f}}$ based on only the fluid-phase-averaged velocity $\boldsymbol{u}_{\mathit{f}}$ , with $\eta _{\mathit{f}}$ the clear-fluid viscosity, would have the very same consequences and not alter the quality of any of the results presented below. Furthermore, we neglect the Reynolds stresses $\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}$ and $\boldsymbol{\sigma ^{\mathit{s}}_{\textit{Re}}}$ , and the generalized-drag force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ merely to keep the notion simple, since they are usually either quadratic in the flow velocity gradients or come with zero-order derivatives, and therefore do not affect linear well-posedness. An exception are virtual-mass forces, which are addressed in § 5.1.1.

Using (5.1) to substitute the respective terms in the buoyancy-improved momentum balances (4.13) and (4.14) yields

(5.2) \begin{align} \partial _t\beta _{\mathit{f}}\rho _{\mathit{f}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{u}_{\mathit{f}}&=0,\end{align}

(5.3) \begin{align} \beta _{\mathit{f}}\rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}&=-\beta _{\mathit{f}}\left (\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}+\mathcal{I}\beta _{\mathit{f}}\underline {\,\,}\right )^{-1}(\boldsymbol{\nabla }P_{\mathit{f}}-\eta _{\mathit{m}}\Delta \boldsymbol{u}_{\mathit{m}})+\beta _{\mathit{f}}\rho _{\mathit{f}}\boldsymbol{g},\end{align}

(5.4) \begin{align} \partial _t\beta _{\mathit{s}}\rho _{\mathit{s}}+\boldsymbol{\nabla }\boldsymbol{\cdot }\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{u}_{\mathit{s}}&=0,\end{align}

(5.5) \begin{align} \beta _{\mathit{s}}\rho _{\mathit{s}}{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}+P_{\mathit{s}}^\prime \boldsymbol{\nabla }\beta _{\mathit{s}}&=-\mathcal{L}\beta _{\mathit{s}}\left (\mathcal{L}\beta _{\mathit{s}}\underline {\,\,}+\mathcal{I}\beta _{\mathit{f}}\underline {\,\,}\right )^{-1}(\boldsymbol{\nabla }P_{\mathit{f}}-\eta _{\mathit{m}}\Delta \boldsymbol{u}_{\mathit{m}})+\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{g},\end{align}

where (5.2) and (5.4) are the fluid and dispersed-phase mass balances (Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026). The model is closed by the volume fraction condition $\beta _{\mathit{f}}+\beta _{\mathit{s}}=1$ and two equations of state, $\rho _{\mathit{s}}=\mathrm{const}$ and $\rho _{\mathit{f}}=\rho _{\mathit{f}}[P_{\mathit{f}}]$ , with $\rho _{\mathit{f}}^\prime \geqslant 0$ .

The above model stands as a simple example for the (typically more complex) two-fluid models used in the chemical engineering community (for an overview, see Lhuillier et al. (Reference Lhuillier, Chang and Theofanous2013)). It is chosen for demonstration purposes, not for physical realism. In fact, if this model is linearly well-posed, so will be most such models as they usually incorporate more dissipative mechanisms. (Non-dissipative virtual-mass forces are addressed in § 5.1.1.) In the analysis below, we will also take a look at the special cases where $\eta _{\mathit{m}}=0$ and/or $P_{\mathit{s}}=0$ . If both are zero, the model resembles the one-pressure model by Stewart & Wendroff (Reference Stewart and Wendroff1984), except for our different buoyancy closure. To transform our model into the latter, one would then only need to apply the small-particle approximation ( $\mathcal{L}\mapsto \mathcal{I}\kern2pt$ ).

To test for linear well-posedness, we consider small-amplitude wave perturbations (subscript ‘ $\ast$ ’) of the dependent variables in (5.2)–(5.5) around a base state (subscript ‘ $0$ ’):

(5.6) \begin{align} &\beta _{\mathit{s}}=\beta _{\mathit{s}0}+\beta _{\mathit{s}\ast },&&\boldsymbol{u}_{\mathit{s}}=\boldsymbol{u}_{\mathit{s}0}+\boldsymbol{u}_{\mathit{s}\ast },&&\boldsymbol{u}_{\mathit{f}}=\boldsymbol{u}_{\mathit{f}0}+\boldsymbol{u}_{\mathit{f}\ast },&&P_{\mathit{f}}=P_{\mathit{f}0}+P_{\mathit{f}\ast }. \end{align}

In the limit of infinitely large wavenumbers and frequencies of these perturbations, relevant for a Hadamard instability, the base state, which varies over finite characteristic length and time, can be assumed to be locally constant (Joseph & Saut Reference Joseph and Saut1990). Then, linearizing (5.2)–(5.5) in these perturbations yields

(5.7) \begin{align} -\rho _{\mathit{f}0}\partial _t\beta _{\mathit{s}\ast }-\rho _{\mathit{f}0}\boldsymbol{u}_{\mathit{f}0}\boldsymbol{\cdot }\boldsymbol{\nabla }\beta _{\mathit{s}\ast }+\beta _{\mathit{f}0}\rho _{\mathit{f}0}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{\mathit{f}\ast }&=-\beta _{\mathit{f}0}\rho _{\mathit{f}0}^\prime \partial _tP_{\mathit{f}\ast }-\beta _{\mathit{f}0}\rho _{\mathit{f}0}^\prime \boldsymbol{u}_{\mathit{f}0}\boldsymbol{\cdot }\boldsymbol{\nabla }P_{\mathit{f}\ast },\end{align}

(5.8) \begin{align} \beta _{\mathit{f}0}\rho _{\mathit{f}0}\partial _t\boldsymbol{u}_{\mathit{f}\ast }+\beta _{\mathit{f}0}\rho _{\mathit{f}0}\boldsymbol{u}_{\mathit{f}0}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{f}\ast }&=-\beta _{\mathit{f}0}\mathcal{N}(\boldsymbol{\nabla }P_{\mathit{f}\ast }-\eta _{\mathit{m}0}\Delta \boldsymbol{u}_{\mathit{m}1\ast })+(\beta _{\mathit{f}}\rho _{\mathit{f}})_1\boldsymbol{g},\end{align}

(5.9) \begin{align} \rho _{\mathit{s}}\partial _t\beta _{\mathit{s}\ast }+\rho _{\mathit{s}}\boldsymbol{u}_{\mathit{s}0}\boldsymbol{\cdot }\boldsymbol{\nabla }\beta _{\mathit{s}\ast }+\beta _{\mathit{s}0}\rho _{\mathit{s}}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}_{\mathit{s}\ast }&=0,\end{align}

(5.10) \begin{align} \beta _{\mathit{s}0}\rho _{\mathit{s}}\partial _t\boldsymbol{u}_{\mathit{s}\ast }+\beta _{\mathit{s}0}\rho _{\mathit{s}}\boldsymbol{u}_{\mathit{s}0}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}_{\mathit{s}\ast }+P_{\mathit{s}0}^\prime \boldsymbol{\nabla }\beta _{\mathit{s}\ast }&=-\beta _{\mathit{s}0}\mathcal{M}^2\mathcal{N}(\boldsymbol{\nabla }P_{\mathit{f}\ast }-\eta _{\mathit{m}0}\Delta \boldsymbol{u}_{\mathit{m}1\ast })+\beta _{\mathit{s}}\rho _{\mathit{s}}\boldsymbol{g},\end{align}

where $\rho _{\mathit{f}0}\equiv \rho _{\mathit{f}}[P_{\mathit{f}0}]$ , $\rho _{\mathit{f}0}^\prime \equiv \rho _{\mathit{f}}^\prime [P_{\mathit{f}0}]$ , $\beta _{\mathit{f}0}\equiv 1-\beta _{\mathit{s}0}$ , $P_{\mathit{s}0}^\prime \equiv P_{\mathit{s}}^\prime [\beta _{\mathit{s}0}]$ , $\eta _{\mathit{m}0}\equiv \eta _{\mathit{m}}[\beta _{\mathit{s}0}]$ , $\mathcal{N}\equiv (\beta _{\mathit{s}0}\mathcal{M}^2+\beta _{\mathit{f}0}\mathcal{I})^{-1}$ , $\boldsymbol{u}_{\mathit{m}1\ast }\equiv \beta _{\mathit{f}0}\boldsymbol{u}_{\mathit{f}\ast }+\beta _{\mathit{s}0}\boldsymbol{u}_{\mathit{s}\ast }-(\boldsymbol{u}_{\mathit{f}0}-\boldsymbol{u}_{\mathit{s}0})\beta _{\mathit{s}\ast }$ and $(\beta _{\mathit{f}}\rho _{\mathit{f}})_1$ denotes the linearization of $\beta _{\mathit{f}}\rho _{\mathit{f}}$ . Note that we used (4.21), assuming identically sized spheres for simplicity (discussed at the end of this section). Now, we Fourier-transform (5.7)–(5.10) in space and time – formally, $\boldsymbol{\nabla }\mapsto \mathrm{i}\boldsymbol{k}$ , $\mathcal{M}\mapsto M_K$ , $\mathcal{N}\mapsto N_K\equiv 1/(\beta _{\mathit{s}0}M_K^2+\beta _{\mathit{f}0})$ and $\partial _t\mapsto -\mathrm{i}\varOmega$ , where $\varOmega$ is the wave frequency – and keep only the terms proportional to $\varOmega$ , $k\equiv |\boldsymbol{k}|$ and $k^2$ ; terms of $O(1)$ , which can potentially also cause Hadamard instabilities, will be considered indirectly via analysis of the dimensionality of the solution space (explained shortly). This formally implies $M_K\mapsto 0$ , $N_K\mapsto 1/\beta _{\mathit{f}0}$ and $\boldsymbol{g}\mapsto 0$ in the Fourier transforms of (5.8) and (5.10). Put together, these steps yield

(5.11) \begin{align} \beta _{\mathit{f}0}\rho _{\mathit{f}0}^\prime (\varOmega -\boldsymbol{u}_{\mathit{f}0}\boldsymbol{\cdot }\boldsymbol{k})\hat P_{\mathit{f}\ast }-\beta _{\mathit{f}0}\rho _{\mathit{f}0}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{\hat u}_{\mathit{f}\ast }-\rho _{\mathit{f}0}(\varOmega -\boldsymbol{u}_{\mathit{f}0}\boldsymbol{\cdot }\boldsymbol{k})\hat \beta _{\mathit{s}\ast }&=0, \\[-12pt]\nonumber \end{align}

(5.12) \begin{align} -\hat P_{\mathit{f}\ast }\boldsymbol{k}+\beta _{\mathit{f}0}\rho _{\mathit{f}0}(\varOmega -\boldsymbol{u}_{\mathit{f}0}\boldsymbol{\cdot }\boldsymbol{k})\boldsymbol{\hat u}_{\mathit{f}\ast }+\mathrm{i}\eta _{\mathit{m}0}k^2\boldsymbol{\hat u}_{\mathit{m}1\ast }&=0, \\[-12pt]\nonumber \end{align}

(5.13) \begin{align} \rho _{\mathit{s}}(\varOmega -\boldsymbol{u}_{\mathit{s}0}\boldsymbol{\cdot }\boldsymbol{k})\hat \beta _{\mathit{s}\ast }-\beta _{\mathit{s}0}\rho _{\mathit{s}}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{\hat u}_{\mathit{s}\ast }&=0, \\[-12pt]\nonumber \end{align}

(5.14) \begin{align} -P_{\mathit{s}0}^\prime \hat \beta _{\mathit{s}\ast }\boldsymbol{k}+\beta _{\mathit{s}0}\rho _{\mathit{s}}(\varOmega -\boldsymbol{u}_{\mathit{s}0}\boldsymbol{\cdot }\boldsymbol{k})\boldsymbol{\hat u}_{\mathit{s}\ast }&=0, \end{align}

where the hat indicates the Fourier transformation with respect to $\boldsymbol{x}$ and $t$ . To simplify the problem, we decompose $\boldsymbol{\hat u}_{\mathit{f}\ast }$ and $\boldsymbol{\hat u}_{\mathit{s}\ast }$ into longitudinal perturbations parallel to $\boldsymbol{k}$ (superscript ‘ $\parallel$ ’) and transverse perturbations normal to $\boldsymbol{k}$ (superscript ‘ $\perp$ ’),

(5.15) \begin{align} \boldsymbol{\hat u}^\parallel _{\mathit{f}\ast }&\equiv \boldsymbol{\hat u}_{\mathit{f}\ast }\boldsymbol{\cdot }\boldsymbol{e_k}\boldsymbol{e_k}\equiv \hat u^\parallel _{\mathit{f}\ast }\boldsymbol{e_k},&\boldsymbol{\hat u}^\parallel _{\mathit{s}\ast }&\equiv \boldsymbol{\hat u}_{\mathit{s}\ast }\boldsymbol{\cdot }\boldsymbol{e_k}\boldsymbol{e_k}\equiv \hat u^\parallel _{\mathit{s}\ast }\boldsymbol{e_k}, \\[-12pt]\nonumber \end{align}

(5.16) \begin{align} \boldsymbol{\hat u}^\perp _{\mathit{f}\ast }&\equiv \boldsymbol{\hat u}_{\mathit{f}\ast }-\boldsymbol{\hat u}^\parallel _{\mathit{f}\ast },&\boldsymbol{\hat u}^\perp _{\mathit{s}\ast }&\equiv \boldsymbol{\hat u}_{\mathit{s}\ast }-\boldsymbol{\hat u}^\parallel _{\mathit{s}\ast }, \end{align}

where $\boldsymbol{e_k}=\boldsymbol{k}/k$ is the unit vector in the direction of $\boldsymbol{k}$ . Furthermore, without loss of generality, we use a Cartesian coordinate system with one of its axes aligned with $\boldsymbol{e_k}$ . Defining

(5.17) \begin{align} \boldsymbol{X}\equiv \big (\hat P_{\mathit{f}\ast },\hat u^\parallel _{\mathit{f}\ast },\hat u^\perp _{\mathit{f}1\ast },\hat u^\perp _{\mathit{f}2\ast },\hat \beta _{\mathit{s}\ast },\hat u^\parallel _{\mathit{s}\ast },\hat u^\perp _{\mathit{s}1\ast },\hat u^\perp _{\mathit{s}2\ast }\big)^{\textrm{T}}, \end{align}

where the indices $1$ and $2$ indicate the two orthogonal transverse directions, we can therefore write (5.11)–(5.14) as

(5.18) \begin{align} \boldsymbol{C}\boldsymbol{X}=0. \end{align}

Here $\boldsymbol{C}$ is a $8\times 8$ matrix given by

(5.19) \begin{align} \boldsymbol{C}\equiv \left (\begin{matrix} \boldsymbol{C}^\parallel _{\mathit{f}} & \boldsymbol{0} & {\boldsymbol{C}^{\mathit{s}}_{\mathit{f}}} & \\[3pt] \boldsymbol{0} & \boldsymbol{C}^\perp _{\mathit{f}} & & \\[3pt] \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{C}^\parallel _{\mathit{s}} & \boldsymbol{0} \\[3pt] \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{0} & \boldsymbol{C}^\perp _{\mathit{s}} \end{matrix}\right )\!, \end{align}

with $\boldsymbol{C}^\parallel _{\mathit{f}}$ , $\boldsymbol{C}^\perp _{\mathit{f}}$ , $\boldsymbol{C}^\parallel _{\mathit{s}}$ and $\boldsymbol{C}^\perp _{\mathit{s}}$ $2\times 2$ matrices defined, respectively, as

(5.20) \begin{align} \boldsymbol{C}^\parallel _{\mathit{f}}&\equiv \begin{pmatrix} \beta _{\mathit{f}0}\rho _{\mathit{f}0}^\prime \big (\varOmega -u^\parallel _{\mathit{f}0}k\big ) & -\beta _{\mathit{f}0}\rho _{\mathit{f}0}k \\[7pt] -k & \beta _{\mathit{f}0}\rho _{\mathit{f}0}\big (\varOmega -u^\parallel _{\mathit{f}0}k+\mathrm{i}\eta _{\mathit{m}}k^2/\rho _{\mathit{f}0}\big ) \end{pmatrix}\!, \\[-12pt]\nonumber \end{align}

(5.21) \begin{align} \boldsymbol{C}^\perp _{\mathit{f}}&\equiv \begin{pmatrix} \beta _{\mathit{f}0}\rho _{\mathit{f}0}\big(\varOmega -u^\parallel _{\mathit{f}0}k+\mathrm{i}\eta _{\mathit{m}}k^2/\rho _{\mathit{f}0}\big ) & 0 \\[7pt] 0 & \beta _{\mathit{f}0}\rho _{\mathit{f}0}\big(\varOmega -u^\parallel _{\mathit{f}0}k+\mathrm{i}\eta _{\mathit{m}}k^2/\rho _{\mathit{f}0}\big ) \end{pmatrix}\!, \\[-12pt]\nonumber \end{align}

(5.22) \begin{align} \boldsymbol{C}^\parallel _{\mathit{s}}&\equiv \begin{pmatrix} \rho _{\mathit{s}}\big (\varOmega -u^\parallel _{\mathit{s}0}k\big ) & -\beta _{\mathit{s}0}\rho _{\mathit{s}}k \\[7pt] -P_{\mathit{s}0}^\prime k & \beta _{\mathit{s}0}\rho _{\mathit{s}}\big (\varOmega -u^\parallel _{\mathit{s}0}k\big ) \end{pmatrix}\!, \\[-12pt]\nonumber \end{align}

(5.23) \begin{align} \boldsymbol{C}^\perp _{\mathit{s}}&\equiv \begin{pmatrix} \beta _{\mathit{s}0}\rho _{\mathit{s}}\big (\varOmega -u^\parallel _{\mathit{s}0}k\big ) & 0 \\[7pt] 0 & \beta _{\mathit{s}0}\rho _{\mathit{s}}\big (\varOmega -u^\parallel _{\mathit{s}0}k\big ) \end{pmatrix}\!, \end{align}

and $\boldsymbol{C}^{\mathit{s}}_{\mathit{f}}$ a $4\times 4$ matrix defined as

(5.24) \begin{align} \boldsymbol{C}^{\mathit{s}}_{\mathit{f}}\equiv \begin{pmatrix} -\rho _{\mathit{f}0}\big (\varOmega -u^\parallel _{\mathit{f}0}k\big ) & 0 & 0 & 0 \\[7pt] -\mathrm{i}\eta _{\mathit{m}0}k^2\big (u^\parallel _{\mathit{f}0}-u^\parallel _{\mathit{s}0}\big ) & \mathrm{i}\beta _{\mathit{s}0}\eta _{\mathit{m}0}k^2 & 0 & 0 \\[7pt] -\mathrm{i}\eta _{\mathit{m}0}k^2\big (u^\perp _{\mathit{f}10}-u^\perp _{\mathit{s}10}\big ) & 0 & \mathrm{i}\beta _{\mathit{s}0}\eta _{\mathit{m}0}k^2 & 0 \\[7pt] -\mathrm{i}\eta _{\mathit{m}0}k^2 \big(u^\perp _{\mathit{f}20}-u^\perp _{\mathit{s}20}\big ) & 0 & 0 & \mathrm{i}\beta _{\mathit{s}0}\eta _{\mathit{m}0}k^2 \end{pmatrix}\!, \end{align}

where the subscript ‘ $0$ ’ indicates again the base state value of a quantity; for example, $u^\perp _{\mathit{f}10}$ is the base state value of $u^\perp _{\mathit{f}1}$ , the first transverse component of $\boldsymbol{u}_{\mathit{f}}$ .

We now look for those frequencies $\varOmega$ that allow for non-trivial perturbation vectors $\boldsymbol{X}\ne 0$ in (5.18), which can be obtained from the condition

(5.25) \begin{align} \mathrm{det}(\boldsymbol{C})=\mathrm{det}(\boldsymbol{C}^\parallel _{\mathit{f}})\mathrm{det}(\boldsymbol{C}^\perp _{\mathit{f}})\mathrm{det}(\boldsymbol{C}^\parallel _{\mathit{s}})\mathrm{det}(\boldsymbol{C}^\perp _{\mathit{s}})=0. \end{align}

It is equivalent to the condition

(5.26) \begin{align} \mathrm{det}(\boldsymbol{C}^\parallel _{\mathit{f}})=0\lor \mathrm{det}(\boldsymbol{C}^\perp _{\mathit{f}})=0\lor \mathrm{det}(\boldsymbol{C}^\parallel _{\mathit{s}})=0\lor \mathrm{det}(\boldsymbol{C}^\perp _{\mathit{s}})=0 \end{align}

and yields eight corresponding frequency solutions:

(5.27) \begin{align} \eta _{\mathit{m}0}&=0:&\varOmega ^\parallel _{\mathit{f}1[2]}&=\big (u^\parallel _{\mathit{f}0}+[-]\big (\beta _{\mathit{f}0}\rho _{\mathit{f}0}^\prime \big )^{-1/2}\big )k,&\varOmega ^\perp _{\mathit{f}1,2}&=u^\parallel _{\mathit{f}0}k, \\[-12pt]\nonumber \end{align}

(5.28) \begin{align} \eta _{\mathit{m}0}&\gt 0:&\varOmega ^\parallel _{\mathit{f}1}&\simeq u^\parallel _{\mathit{f}0}k,&\varOmega ^\parallel _{\mathit{f}2},\varOmega ^\perp _{\mathit{f}1,2}&\simeq -\mathrm{i}\eta _{\mathit{m}0}k^2/\rho _{\mathit{f}0}, \\[-12pt]\nonumber \end{align}

(5.29) \begin{align} \eta _{\mathit{m}0}&\geqslant 0:&\varOmega ^\parallel _{\mathit{s}1[2]}&=\big (u^\parallel _{\mathit{s}0}+[-]\big (P_{\mathit{s}0}^\prime /\rho _{\mathit{s}}\big )^{1/2}\big )k,&\varOmega ^\perp _{\mathit{s}1,2}&=u^\parallel _{\mathit{s}0}k. \end{align}

Here the symbol ‘ $\simeq$ ’ denotes the asymptotic expansion in the limit $k\to \infty$ , the subscript ‘ $1,2$ ’ indicates a double solution, while the subscript ‘ $1[2]$ ’ means that the second solution refers to the symbol surrounded by square brackets on the corresponding right-hand side, whereas the first solution refers to the symbol immediately before the square brackets. Note that, for an incompressible fluid, $\rho ^\prime _{\mathit{f}0}=0$ , the determinant $\mathrm{det}(\boldsymbol{C}^\parallel _{\mathit{f}})=-\beta _{\mathit{f}0}\rho _{\mathit{f}0}k^2$ is independent of $\varOmega$ and strictly smaller than zero, implying that the double solution $\varOmega ^\parallel _{\mathit{f}1,2}$ no longer exist, while the other solutions remain unaffected.

Due to $\rho _{\mathit{f}0}^\prime \geqslant 0$ and $P^\prime _{\mathit{s}}\geqslant 0$ , all the above frequency solutions are asymptotically either negative imaginary, associated with damping or real-valued, associated with pure oscillations, implying that the $O(k)$ and $O(k^2)$ terms themselves are unable to cause a Hadamard instability. However, the system of PDEs is not necessarily linearly well-posed, since the previously neglected $O(1)$ terms can conspire with the $O(k)$ terms to cause unbounded growth as $k^b$ , with $0\lt b\lt 1$ (Langham et al. Reference Langham, Meng, Webb, Johnson and Gray2025b ). This happens when $n\geqslant 2$ real-valued frequency solutions asymptotically coincide for a given base state and their corresponding null spaces (i.e. the spaces of the non-trivial perturbation vector solutions) are less than $n$ -dimensional (Langham et al. Reference Langham, Meng, Webb, Johnson and Gray2025b ). To check for this possibility we calculate the basis vectors of the null spaces corresponding to the frequency solutions in (5.27)–(5.29):

(5.30) \begin{align} \eta _{\mathit{m}0}&=0:&\boldsymbol{X}^\parallel _{\mathit{f}1[2]}&=\big (+[-]\rho _{\mathit{f}0}\big (\beta _{\mathit{f}0}/\rho _{\mathit{f}0}^\prime \big )^{1/2},1,0,0,0,0,0,0\big )^{\textrm{T}}, \\[-12pt]\nonumber \end{align}

(5.31) \begin{align} \eta _{\mathit{m}0}&\gt 0:&\boldsymbol{X}^\parallel _{\mathit{f}1[2]}&\simeq \left (1[0],0[1],0,0,0,0,0,0\right )^{\textrm{T}}, \\[-12pt]\nonumber \end{align}

(5.32) \begin{align} \eta _{\mathit{m}0}&\geqslant 0:&\boldsymbol{X}^\perp _{\mathit{f}1[2]}&=\left (0,0,1[0],0[1],0,0,0,0\right )^{\textrm{T}}, \\[-12pt]\nonumber \end{align}

(5.33) \begin{align} \eta _{\mathit{m}0}&=0:&\boldsymbol{X}^\parallel _{\mathit{s}1[2]}&=\left (\beta _{\mathit{f}0}\rho _{\mathit{f}0}\big (u^\parallel _{\mathit{s}0}+[-]\big (P_{\mathit{s}0}^\prime /\rho _{\mathit{s}}\big )^{1/2}-u^\parallel _{\mathit{f}0}\big ),1,0, \right. \nonumber \\ &&&\quad \;\; \left. 0,\chi _{1[2]},+[-]\beta _{\mathit{s}0}^{-1}\big (P_{\mathit{s}0}^\prime /\rho _{\mathit{s}}\big )^{1/2}\chi _{1[2]},0,0\right )^{\textrm{T}}, \\[-12pt]\nonumber \end{align}

(5.34) \begin{align} \eta _{\mathit{m}0}&\gt 0:&\boldsymbol{X}^\parallel _{\mathit{s}1[2]}&\simeq \left (0,\beta _{\mathit{s}0}\left (u^\parallel _{\mathit{f}0}-u^\parallel _{\mathit{s}0}-[+]\left (P_{\mathit{s}0}^\prime /\rho _{\mathit{s}}\right )^{1/2}\right ),\beta _{\mathit{s}0}\big (u^\perp _{\mathit{f}10}-u^\perp _{\mathit{s}10}\big ),\right . \nonumber \\ &&&\quad \;\;\left .\beta _{\mathit{s}0}\big (u^\perp _{\mathit{f}20}-u^\perp _{\mathit{s}20}\big ),\beta _{\mathit{f}0}\beta _{\mathit{s}0},+[-]\beta _{\mathit{f}0}\left (P_{\mathit{s}0}^\prime /\rho _{\mathit{s}}\right )^{1/2},0,0\right )^{\textrm{T}}\!, \\[-12pt]\nonumber \end{align}

(5.35) \begin{align} \eta _{\mathit{m}0}&=0:&\boldsymbol{X}^\perp _{\mathit{s}1[2]}&=(0,0,0,0,0,0,1[0],0[1])^{\textrm{T}}, \\[-12pt]\nonumber \end{align}

(5.36) \begin{align} \eta _{\mathit{m}0}&\gt 0:&\boldsymbol{X}^\perp _{\mathit{s}1[2]}&\simeq (0,0,-\beta _{\mathit{s}0}/\beta _{\mathit{f}0}[0],0[-\beta _{\mathit{s}0}/\beta _{\mathit{f}0}],0,0,1[0],0[1])^{\textrm{T}}, \end{align}

where

(5.37) \begin{align} \chi _{1[2]}\equiv -\beta _{\mathit{f}0}\left (\Big (u^\parallel _{\mathit{s}0}+[-]\big (P_{\mathit{s}0}^\prime /\rho _{\mathit{s}}\big )^{1/2}-u^\parallel _{\mathit{f}0}\Big )^{-1}-\beta _{\mathit{f}0}\rho _{\mathit{f}0}^\prime \Big (u^\parallel _{\mathit{s}0}+[-]\big (P_{\mathit{s}0}^\prime /\rho _{\mathit{s}}\big )^{1/2}-u^\parallel _{\mathit{f}0}\Big )\right )\!. \end{align}

One immediately obtains that the basis vectors corresponding to the transverse modes, $\boldsymbol{X}^\perp _{\mathit{f}1}$ , $\boldsymbol{X}^\perp _{\mathit{f}2}$ , $\boldsymbol{X}^\perp _{\mathit{s}1}$ and $\boldsymbol{X}^\perp _{\mathit{s}2}$ , are always linearly independent from each other and from the basis vectors corresponding to the longitudinal modes, i.e. they never cause problems. Furthermore, we can infer the following statements.

1. (i) If variations of the effective dispersed-phase pressure are neglected, $P^\prime _{\mathit{s}}=0$ , then $\varOmega ^\parallel _{\mathit{s}1}=\varOmega ^\parallel _{\mathit{s}2}$ and $\boldsymbol{X}^\parallel _{\mathit{s}1}=\boldsymbol{X}^\parallel _{\mathit{s}2}$ regardless of the behaviour of $\eta _{\mathit{m}}$ . Hence, the null space corresponding to $\varOmega ^\parallel _{\mathit{s}1}$ and $\varOmega ^\parallel _{\mathit{s}2}$ is just one-dimensional rather than two-dimensional, implying that the system of PDEs is ill-posed due to unbounded growth as $k^{1/2}$ (Langham et al. Reference Langham, Meng, Webb, Johnson and Gray2025b ).

2. (ii) If $P_{\mathit{s}}^\prime \gt 0$ and if the fluid is incompressible, $\rho _{\mathit{f}0}^\prime =0$ , then the system is linearly well-posed regardless of the behaviour of $\eta _{\mathit{m}}$ , since the frequency solutions $\varOmega ^\parallel _{\mathit{f}1}$ and $\varOmega ^\parallel _{\mathit{f}2}$ and their corresponding basis vectors $\boldsymbol{X}^\parallel _{\mathit{f}1}$ and $\boldsymbol{X}^\parallel _{\mathit{f}2}$ no longer exist, while the remaining six basis vectors are always linearly independent from each other.

3. (iii) If $P_{\mathit{s}}^\prime \gt 0$ and if the fluid is compressible, $\rho _{\mathit{f}0}^\prime \gt 0$ , then the system is linearly well-posed for $\eta _{\mathit{m}}\gt 0$ but generally ill-posed if viscosity effects are neglected, $\eta _{\mathit{m}}=0$ . In the latter case, for those base states that satisfy $\varOmega ^\parallel _{\mathit{f}1}=\varOmega ^\parallel _{\mathit{s}1}$ or $\varOmega ^\parallel _{\mathit{f}1}=\varOmega ^\parallel _{\mathit{s}2}$ or $\varOmega ^\parallel _{\mathit{f}2}=\varOmega ^\parallel _{\mathit{s}1}$ or $\varOmega ^\parallel _{\mathit{f}2}=\varOmega ^\parallel _{\mathit{s}2}$ , the corresponding basis vectors coincide, i.e. $\boldsymbol{X}^\parallel _{\mathit{f}1}=\boldsymbol{X}^\parallel _{\mathit{s}1}$ or $\boldsymbol{X}^\parallel _{\mathit{f}1}=\boldsymbol{X}^\parallel _{\mathit{s}2}$ or $\boldsymbol{X}^\parallel _{\mathit{f}2}=\boldsymbol{X}^\parallel _{\mathit{s}1}$ or $\boldsymbol{X}^\parallel _{\mathit{f}2}=\boldsymbol{X}^\parallel _{\mathit{s}2}$ , respectively, implying again unbounded growth as $k^{1/2}$ (Langham et al. Reference Langham, Meng, Webb, Johnson and Gray2025b ). By contrast, viscous dissipation for $\eta _{\mathit{m}}\gt 0$ prevents such resonances from occurring through asymptotically decoupling fluid phase and dispersed-phase oscillations from one another.

Note that, if we modelled viscous dissipation as $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau ^{\,\mathit{f}}}=\eta _{\mathit{f}}\Delta \boldsymbol{u}_{\mathit{f}}$ instead of $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau ^{\,\mathit{f}}}=\eta _{\mathit{m}}\Delta \boldsymbol{u}_{\mathit{m}}$ , the above statements would remain true if $\eta _{\mathit{m}}$ were formally replaced by $\eta _{\mathit{f}}$ .

In summary, even simplistic two-fluid models are linearly well-posed when employing our buoyancy closure. One only requires standard viscous damping of the fluid phase and an effective dispersed-phase pressure $P_{\mathit{s}}$ that strictly monotonously increases with the dispersed-phase volume fraction $\beta _{\mathit{s}}$ , which is a requirement dictated by physics (e.g. granular kinetic theory (Berzi Reference Berzi2024)).

In the above analysis, we have argued on the basis of (4.21), which is valid only for systems of identically sized spheres. However, it is straightforward to generalize this analysis to arbitrary particle sizes and shapes, giving the same results. This is because the very definition of the operator $\mathcal{L}$ in (4.1) always means that spatial variations smaller than the characteristic particle size are smoothed out, attenuating buoyancy effects as $k^{-4}$ when $k\to \infty$ .

5.1.1. Virtual-mass force density

In addition to buoyancy, virtual-mass-related effects have also been linked to the ill-posedness of two-fluid models (Lhuillier et al. Reference Lhuillier, Chang and Theofanous2013). A standard expression for the virtual-mass force density $\beta _{\mathit{s}}\langle \boldsymbol{f}_{\mathit{v}}\rangle ^{\mathit{s}}\equiv \beta _{\mathit{s}}\langle \sum _pX^p\boldsymbol{F}^p_{\mathit{v}}/V^p\rangle ^{\mathit{s}}$ , where $\boldsymbol{F}^p_{\mathit{v}}$ is the virtual-mass force acting on a particle $p$ , is that proposed by Jackson (Reference Jackson2000):

(5.38) \begin{align} \beta _{\mathit{s}}\left \langle \boldsymbol{f}_{\mathit{v}}\right \rangle ^{\mathit{s}}=C_{\mathit{v}}\rho _{\mathit{f}}\beta _{\mathit{s}}\beta _{\mathit{f}}\big ({\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}\big ), \end{align}

where $C_{\mathit{v}}[\beta _{\mathit{s}}]$ is the virtual-mass coefficient, with $C_{\mathit{v}}[0]=1/2$ . We have written (5.38) in a form that involves multiplication with $\beta _{\mathit{f}}$ , analogous to the generalized-drag force density expression in § 3.2, in order to highlight that $\beta _{\mathit{s}}\left \langle \boldsymbol{f}_{\mathit{v}}\right \rangle ^{\mathit{s}}$ is a force per unit volume of the fluid–particle mixture. Incorporating (5.38) as an interaction force density in (5.3) and (5.5) results in frequency solutions $\varOmega \sim k$ with positive imaginary part for many base states and therefore in growth as $k^1$ and ill-posedness. However, if one took into account that $\boldsymbol{u}_{\mathit{f}}$ , which plays a role analogous to the background flow velocity $\boldsymbol{\tilde u}$ in the MRG equation (1.1), should actually be averaged over the particles’ volumes (e.g. see (1.1)), one should model $\beta _{\mathit{s}}\left \langle \boldsymbol{f}_{\mathit{v}}\right \rangle ^{\mathit{s}}$ as

(5.39) \begin{align} \beta _{\mathit{s}}\left \langle \boldsymbol{f}_{\mathit{v}}\right \rangle ^{\mathit{s}}=C_{\mathit{v}}\rho _{\mathit{f}}\beta _{\mathit{f}}\big (\mathcal{L}\beta _{\mathit{s}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\beta _{\mathit{s}}{\mathrm{D}}^{\mathit{s}}_t\boldsymbol{u}_{\mathit{s}}\big ). \end{align}

This becomes immediately clear when comparing with the manner in which the average generalized-buoyancy force on a particle $p$ , $\langle \boldsymbol{F}^p_{\mathit{B}}\rangle =\int _{\mathbb{V}^p}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\hat \sigma ^{\,\mathit{f}}}{\mathrm{d}} V$ is linked to the corresponding generalized-buoyancy force density $\beta _{\mathit{s}}\left \langle \boldsymbol{f}_{\!\mathit{B}}\right \rangle ^{\mathit{s}}$ in (4.3). Incorporating (5.39) as an interaction force density in (5.3) and (5.5) does not alter linear well-posedness, since all above frequency solutions and basis vectors remain the same if the dispersed-phase mass density $\rho _{\mathit{s}}$ is formally replaced as

(5.40) \begin{align} \rho _{\mathit{s}}\mapsto \rho _{\mathit{s}}+C_{\mathit{v}0}\rho _{\mathit{f}0}\beta _{\mathit{f}0}, \end{align}

where $C_{\mathit{v}0}\equiv C_{\mathit{v}}[\beta _{\mathit{s}0}]$ .

5.2. Numerical implementation in two-fluid models

The relevance of our improved buoyancy closure for numerical two-fluid models depends on the size of the mesh grid cells discretizing a given system. Since a value of a physical quantity assigned to a grid cell already represents an average over the cell domain, the effect of the smoothing operator $\mathcal{L}$ distinguishing our closure from its small-particle approximation, given by the closure by Revil-Baudard & Chauchat (Reference Revil-Baudard and Chauchat2013), (2.19), is negligible if the grid size is much larger than $R$ . Hence, for simulations of large-scale systems with two fluid models, it suffices to employ this small-particle approximation. However, if the grid size is of the order of $R$ or smaller, the effects of $\mathcal{L}$ become relevant. Some researchers are under the impression that it makes no physical sense to analyse two-fluid continuum models at such small scales. For example, Jackson (Reference Jackson2000, pp. 112–115), when analysing the stability of homogeneous fluidized beds, stated that, for wave disturbances of the macroscopic flow quantities with wavelengths of the order of $R$ , ‘continuum equations are clearly inappropriate’. However, we disagree with this perception. When deriving two-fluid continuum equations via ensemble averaging the microscopic equations of motion of the fluid and particles, no separation between the sizes of the microscopic and macroscopic scales is mathematically required (Prosperetti Reference Prosperetti1998; Pähtz et al. Reference Pähtz, Chen, Zhu, Tholen and He2026). More importantly, there are, indeed, disperse two-phase flows where macroscopic quantities change over distances (much) smaller than $R$ . A notable example is fluid-driven sediment transport, which is characterized by a typically sudden transition from a densely packed bed to a potentially very dilute transport layer (Chassagne et al. Reference Chassagne, Bonamy and Chauchat2023; Tholen et al. Reference Tholen, Pähtz, Kamath, Parteli and Kroy2023; Fry et al. Reference Fry, Lacaze, Bonometti, Elyakime and Charru2024), e.g. see figures 2 and 3.

When numerical mesh grid cells are of the order of $R$ or smaller, we suggest to use the form of the two-fluid momentum balances in (4.15) and (4.16) to implement our buoyancy closure in two-fluid models. This form involves no modification of the effective fluid phase stress tensor $\boldsymbol{\sigma ^{\,\mathit{f}}}$ , while the operator $\mathcal{L}$ is applied on only the single term $\beta _{\mathit{s}}\boldsymbol{A}$ , with $\boldsymbol{A}\equiv \rho _{\mathit{f}}{\mathrm{D}}^{\,\mathit{f}}_t\boldsymbol{u}_{\mathit{f}}-\rho _{\mathit{f}}\boldsymbol{g}-\beta _{\mathit{f}}^{-1}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\sigma ^{\,\mathit{f}}_{\textit{Re}}}+\beta _{\mathit{f}}^{-1}\beta _{\mathit{s}}\langle \boldsymbol{f}_{\!\mathit{D}}\rangle ^{\mathit{s}}$ . Then, despite the mathematical complexity involved in $\mathcal{L}$ , there is a simple approximate manner with which its effects can be accounted for,

(5.41) \begin{align} \mathcal{L}\beta _{\mathit{s}}\boldsymbol{A}\approx \beta _{\mathit{s}}\mathcal{M}^2\boldsymbol{A}, \end{align}

where the radius $R$ involved in the definition of $\mathcal{M}$ should, in the case of arbitrarily sized and shaped particles, be interpreted as the volume-equivalent radius corresponding to the average particle volume. The approximation (5.41) preserves the regularization effect of $\mathcal{L}$ and, in the case of identically sized spheres, agrees with the exact expression for both extreme limits in which either of the two quantities $\beta _{\mathit{s}}$ and $\boldsymbol{A}$ spatially varies much more than the other one. The numerical implementation of $\mathcal{M}^2\boldsymbol{A}$ is simple. Presume $\boldsymbol{A}$ is given for every numerical grid cell $\boldsymbol{x}_i$ . Then, assign to each grid cell the value of $\boldsymbol{A}$ averaged over a ball of radius $R$ centred at $\boldsymbol{x}_i$ , giving $\mathcal{M}\boldsymbol{A}[\boldsymbol{x}_i]$ . Repeating the procedure for $\mathcal{M}\boldsymbol{A}[\boldsymbol{x}_i]$ results in $\mathcal{M}^2\boldsymbol{A}[\boldsymbol{x}_i]$ .