1. Introduction
The stability of a quiescent fluid layer under the influence of a settling dispersed phase composed of particles, drops or bubbles is a fluid dynamic problem of remarkable richness and complexity (Guazzelli & Hinch Reference Guazzelli and Hinch2011; Mudde Reference Mudde2005). Even under highly idealised physical conditions, the problem statement involves a large number of physical parameters required to specify the material properties of both the fluid and the suspension, as well as all their boundary conditions. This is even more true when the couplings between the fluid and particles are not only mechanical but also energetic, involving, for example, temperature and phase changes (melting, condensation, evaporation). In the case where the particles are small, very numerous, and highly diluted, an Eulerian two-fluid approach can be adopted to identify the parametric conditions under which they can destabilise the fluid, leading to large-scale advective motion (Saffman Reference Saffman1962).
Recently, Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021) proposed studying a system in which a fluid layer is confined between two horizontal plates maintained at different temperatures, using heating from below to create an unstable density stratification (Rayleigh–Bénard configuration). In this system, particles heavier than the fluid are introduced continuously from the top wall, at their terminal velocity and prescribed temperature, and removed from the bottom wall. This system, named particulate Rayleigh–Bénard (pRB), is shown to be more stable than the particle-free system. In other words, the Rayleigh number of the system must be larger than
$Ra_c \simeq 1708$
to trigger large-scale fluid motion, where
$Ra_c$
identifies the supercritical bifurcation point of a pure fluid system. Subsequently, Raza, Hirata & Calzavarini (Reference Raza, Hirata and Calzavarini2024) extended the pRB model to particles of arbitrary density, including particles lighter than the fluid, which are injected from the lower plate at fixed temperatures. Even in this case, the system with particles is more stable than the particle-free system, regardless of the fluid-to-particle relative mass density and the strength of the mechanical and/or thermal couplings between the dispersed and continuous phases. The latter result appears to contradict the analysis of Nakamura et al. (Reference Nakamura, Yoshikawa, Tasaka and Murai2020, Reference Nakamura, Yoshikawa, Tasaka and Murai2021), who demonstrated that the injection of bubbles from the bottom into an isothermal fluid layer is linearly unstable. Differently from the previously mentioned studies, these authors had considered the possibility of injecting bubbles at sub-terminal velocities, which is more realistic with respect to experiments. This fact has an important technical consequence in the linear stability problem – it allows the particle concentration to develop spatial inhomogeneities, which are prevented when the particles are injected at terminal velocity. Nakamura et al. (Reference Nakamura, Yoshikawa, Tasaka and Murai2020) observed that the variation of the injection velocity does not affect the system stability, unless the velocity is very close to the terminal rising velocity.
These issues raise the question of the role of particle injection velocity in the pRB system, and whether there exists specific combinations of injection velocities and particle densities allowing the system to be destabilised or stabilised and, thus, controlled. They are addressed in the present study by improving the linear stability analysis of Raza et al. (Reference Raza, Hirata and Calzavarini2024), where the particle injection velocity was set equal to the terminal velocity. Doing so shows that the system can be either stabilised or destabilised by properly tuning the particle inlet velocity and flux, where heavy and light particles induce opposite trends.
2. The pRB model system
Following Raza et al. (Reference Raza, Hirata and Calzavarini2024), we adopt an Eulerian model to describe the dynamics of a macroscopic material particle suspension in the pRB setting. The particle volume concentration is assumed to be small throughout, allowing the fluid to be treated as incompressible and governed by the Navier–Stokes equations under the Boussinesq approximation for the velocity
$\boldsymbol {u}(\boldsymbol{x},t)$
and temperature
$T(\boldsymbol{x},t)$
fields. However, due to the total conservation of momentum and thermal energy, the particulate phase exerts both mechanical and thermal feedback on the fluid. This phase is characterised by the individual material properties of the particles, including mass density
$\rho _p$
, diameter
$d_p$
and specific heat capacity at constant pressure
$c_{Pp}$
. Additionally, it is described by the volume concentration
$\alpha (\boldsymbol{x},t)$
, velocity
$\boldsymbol {w}(\boldsymbol{x},t)$
and temperature
$T_p(\boldsymbol{x},t)$
fields. The governing conservation equations for mass, momentum and heat for both the fluid and particle phases are given as follows:
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\,\,\, 0 = \boldsymbol {\nabla \cdot u}, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\quad d_t \alpha = - \alpha (\boldsymbol {\nabla \cdot w}), \end{align}
\begin{align}& D_t\boldsymbol{u} = \dfrac {-\nabla p}{\rho } +\nu \nabla ^2 \boldsymbol{u} + [1-\beta _T(T-T_r)]\boldsymbol{g} + \alpha \bigg [ \left (D_t\boldsymbol{u} - \boldsymbol{g}\right ) + \frac {\rho _p}{\rho } \left (\boldsymbol{g} - d_t\boldsymbol{w} \right ) \bigg], \end{align}
\begin{align}&\qquad\qquad\qquad\quad\quad {d}_t \boldsymbol{w} = \beta D_t\boldsymbol{u}+\tau _p^{-1}(\boldsymbol{u}-\boldsymbol{w})+(1-\beta )\boldsymbol{g}, \end{align}
\begin{align}&\qquad\qquad\qquad\quad D_t T = \kappa \nabla ^2 T + \alpha \big[ {D_t T} - E\ \tau _{T}^{-1}(T - T_p)\big], \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad {d}_t T_p= \tau _{T}^{-1}(T-T_p). \end{align}
A few additional observations are in order. First,
$ D_t() = \partial _t () + \boldsymbol {u \cdot \nabla} ()$
and
$ d_t() = \partial _t() + \boldsymbol {w \cdot \nabla} ()$
denote, respectively, the fluid-phase and particulate-phase material derivative operators, where
$\boldsymbol{\nabla} ()$
is the operator containing the spatial derivative vector. Second, we account for three primary hydrodynamic forces acting on the particles: (i) the Stokes drag force, (ii) the fluid acceleration force with the added mass correction and (iii) buoyancy. The drag force is parameterised by the viscous response time
$ \tau _p = d_p^2/(12 \nu \beta ),$
where
$\nu$
is the fluid kinematic viscosity. The intensity of the added mass force is modulated by the modified density ratio
$ \beta = 3\rho /(\rho + 2\rho _p),$
where
$\rho$
denotes the fluid mass density. Third, we assume that the temperature within each particle remains uniform (lumped approximation), with its relaxation towards equilibrium characterised by the time scale
$ \tau _T = d_p^2 E/(12 \kappa ),$
where
$\kappa$
is the fluid thermal diffusivity and
$ E = \rho _p c_{Pp}/(\rho c_P),$
with
$c_P$
representing the fluid specific heat capacity at constant pressure. Finally, the remaining constants include the fluid volumetric thermal expansion coefficient
$\beta _T$
at the reference temperature
$T_r$
, the gravitational acceleration vector
$\textbf{g} = - g \hat{\textbf{z}}$
(with
$\textbf{z}$
vertical unit vector), and the pressure field
$p(\boldsymbol {x},t)$
.
The domain considered here is three-dimensional, infinitely long in both horizontal coordinates which are vertically bounded by two parallel horizontal walls, with distance
$H$
, located at
$ z = \pm H/2$
, with the unit vector
$ \hat {\boldsymbol {z}}$
pointing upwards. The fluid satisfies no-slip boundary conditions at both walls (
$\boldsymbol {u} = 0$
), which are kept isothermal, with a thermal gap
$\Delta T$
. The bottom wall is the warmest, creating an unstable density stratification when
$ \beta _T \gt 0$
. Particles are introduced from one of the boundaries at a constant volumetric flux and with a prescribed velocity
$ \boldsymbol{w}^*$
, expressed as a multiple of the reference terminal velocity, i.e.
$ \boldsymbol{w}_T = (1-\beta )\tau _p \boldsymbol{g}.$
Particles with
$ \beta \lt 1$
, hereafter referred to as ‘heavy particles’, are injected from the top boundary, whereas those with
$ \beta \gt 1$
, i.e. ‘light particles’, are injected from the bottom boundary. The particle inlet temperature is also set to a fixed value,
$ T_p^*$
, to be specified later. Particle accumulation at the opposite boundary is neglected, i.e. they are assumed to be removed from the domain as soon as they reach the opposite wall. It is worth noting that the equations governing the particulate phase are of first order in space, as they lack a dissipation term in the form of a Laplacian. Consequently, only an inlet particulate-phase boundary condition is required to determine their solution.
2.1. Dimensionless system
In preparation for the upcoming stability study, we rewrote the model in dimensionless form using its characteristic height,
$H$
, conductive time scale,
$\mathcal{T} = t \kappa /H^2$
, and fluid density,
$\rho$
. Accordingly, we define the following dimensionless fields:
\begin{equation} \textbf{U}=\boldsymbol{u}\dfrac {H}{\kappa }, \quad {P = \dfrac {p H^2}{\rho \kappa ^2 }}, \quad \Theta = \dfrac {T-T_r}{\Delta T}, \quad \textbf{W}=\boldsymbol{w}\dfrac {H}{\kappa }, \quad \Theta _p = \dfrac {T_{p}-T_r}{\Delta T}, \end{equation}
which represent the dimensionless fluid velocity, pressure and temperature, as well as the particulate velocity and temperature. Maintaining the same notation for the fluid and particle dimensionless material derivatives, i.e.
$D_{\mathcal{T}}$
and
$d_{\mathcal{T}}$
, and differential operator (
$\boldsymbol{\nabla}$
), and explicitly incorporating the particle momentum feedback, the governing equations (2.1)–(2.6) become
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad\quad\,\, 0 = \boldsymbol{\nabla \cdot} \ \textbf{U}, \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\qquad\qquad d_{\mathcal{T}} \alpha = - \alpha (\boldsymbol{\nabla \cdot} \ \textbf{W}), \end{align}
\begin{align}& D_{\mathcal{T}}\textbf{U} = -\nabla P + Pr \nabla ^2 \textbf{U} + Pr Ra \Theta \hat {\textbf{Z}} \nonumber \\ &\qquad\quad\, + \dfrac {\alpha }{2} \big[ (\beta -1) \big( D_{\mathcal{T}}\textbf{U}+ {Ga Pr^{2} \hat {\textbf{Z}}}\big) - 12 Pr (3-\beta )\Phi ^{-2}( \textbf{U} - \textbf{W} ) \big], \end{align}
\begin{align}&\qquad\quad d_{\mathcal{T}} \textbf{W} = \beta \big( D_{\mathcal{T}}\textbf{U} + 12 Pr \ \Phi ^{-2} \left (\textbf{U} - \textbf{W}\right ) \big) - {(1 - \beta )} {Ga Pr^{2} \hat {\textbf{Z}}}, \end{align}
\begin{align}&\qquad\qquad\qquad D_{\mathcal{T}}\Theta = {\nabla ^2 \Theta } + \alpha \big [{D_{\mathcal{T}}\Theta }- 12\Phi ^{-2}\left (\Theta - \Theta _p \right )\big], \end{align}
\begin{align}&\qquad\qquad\qquad\qquad\quad d_{\mathcal{T}} \Theta _p = 12 (E \Phi ^{2})^{-1}\left (\Theta - \Theta _p \right ), \end{align}
with dimensionless boundary conditions
\begin{align}& \textbf{U} = 0,\quad \Theta = 1 \quad \textrm {at} \quad Z=- \tfrac {1}{2}, \quad \textrm {and}\quad \textbf{U} = 0;\quad \Theta = 0 \quad \textrm {at} \quad Z= \tfrac {1}{2}, \end{align}
\begin{align}&\quad \textbf{W} = \textbf{W}^* =W^*\hat {\textbf{Z}},\qquad \alpha = \mathcal{J}/||\textbf{W}^*||, \quad \Theta _p =\Theta _p^* \quad \textrm {at} \quad Z=Z^*, \end{align}
where
$\hat{\textbf{Z}}$
stands for the vertical unit vector and
$Z^*$
denotes the location of the inlet horizontal wall. Note that
$W^*$
can be either positive or negative, depending on whether the particle is lighter or heavier than the fluid, while the inlet flux
$\mathcal{J}$
is defined as always positive. In the equations (2.8)–(2.13) we have introduced the following dimensionless characteristic parameters:
\begin{equation} Ra=\frac {\beta _T \Delta T g H^3}{\nu \kappa }, \quad Pr = \frac {\nu }{\kappa }, \quad Ga = \frac {g H^3}{\nu ^2}, \quad \Phi = \frac {d_p}{H}, \end{equation}
where
$ Ra$
denotes the Rayleigh number, which quantifies the relative strength of thermally induced buoyancy against mechanical and thermal dissipation,
$ Pr$
the Prandtl number, which characterises the fluid-phase diffusive material properties, and
$ Ga$
the Galileo number, which represents the balance between gravitational and viscous forces. Although
$Ga$
does not depend on the particle properties, it becomes a relevant control parameter whenever the coupling between the particle and the fluid is taken into account. In dimensionless units, the terminal velocity expression is
$\textbf{W}_T ={((1-\beta) }/{\beta })\,\,{ (\Phi ^2}/{12})\, Ga Pr \hat {\textbf{Z}}$
. Along with the modified fluid-to-particle density ratio,
$ \beta$
, the volumetric particle inlet flux,
$\mathcal{J}$
, and the inlet velocity and temperatures
$(W^*, \Theta _p^*)$
, they define the full set of control parameters. In total, the model is governed by nine parameters: three associated with the fluid phase
$( Ra, Pr, Ga )$
and six with the particulate phase
$( \Phi , \beta , E, \mathcal{J} , W^*, \Theta _p^* )$
.
3. Linear stability analysis
3.1. Conductive state
In order to identify the onset of natural convection in the pRB system, its equilibrium solution must first be defined. The equilibrium solution chosen here is its steady state with particles settling or rising in a quiescent and conductive fluid; which is expected only to be a function of
$Z$
and that will be denoted with a zero subscript in the following. Hence, one can impose
$\textbf{U}=0$
,
$\textbf{W}=W_0(Z) \hat {\textbf{Z}}$
,
$\alpha =\alpha _0(Z)$
,
$\Theta =\Theta _0(Z)$
and
$\Theta _p=\Theta _{p0}(Z)$
, and rewrite (2.8)–(2.13) as
\begin{align}&\qquad\qquad\qquad\quad D(\alpha _0 W_0) = 0, \end{align}
\begin{align}& W_0 DW_0 = -{12 Pr \beta W_0}{\Phi ^{-2}} - (1-\beta ) {Ga Pr^{2} \hat {\textbf{Z}}}, \end{align}
\begin{align}&\qquad D^2 \Theta _0 - {12 \alpha _0}{\Phi ^{-2}}(\Theta _0-\Theta _{p0}) = 0 , \end{align}
\begin{align}&\quad\,\,\, W_0 D\Theta _{p0} ={12}({E \Phi ^2})^{-1} (\Theta _0-\Theta _{p0}), \end{align}
where
$D$
represents the derivative with respect to
$Z$
. The steady state governing (3.1)–(3.4), with their here respective boundary conditions derived from (2.14) and (2.15), are computed numerically. Particle volume fraction profiles are shown in figure 1. It is possible to note that, in general, particle concentration varies with height, except for the special case where particles are injected at their terminal velocity (
$W^*=W_T$
), which is the case previously studied by Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021) and Raza et al. (Reference Raza, Hirata and Calzavarini2024). When particles are introduced at a velocity lower (higher) than their terminal velocity, they accelerate (decelerate) until they reach it. This process leads to an accumulation (rarefaction) of particles near the injection wall, e.g. in the upper region for heavier particles and in the lower region for lighter ones. Finally, for the prescribed particle inflow considered here, reducing the injection velocity increases the volumetric concentration of incoming particles. Hence, smaller inlet velocities lead to a higher amount of particles.
Vertical dependence of the steady-state volume concentration field of the particulate phase,
$\alpha _0(Z)$
, for different inlet velocities
$W^*$
for (a)
$\beta =0.5$
(heavy particles) and (b)
$\beta =3$
(bubbles). Data are normalised by the steady-state concentration field corresponding to the case
$W^*=W_T$
.
3.2. Perturbation equations
Having defined the equilibrium solution of the pRB system, the next step is to analyse its linear stability. This is done here by decomposing all variables into their respective steady-state and small-amplitude perturbations (denoted with a prime index), linearising the resulting equations and assuming these perturbations be modelled as normal modes, such as
\begin{equation} \xi '(X,Z,\mathcal{T})=\xi ^n(Z) \exp (i k X+\lambda \mathcal{T}) + \textit {c.c.} \end{equation}
where
$\xi '=\{\textbf{U}',\textbf{W}', \alpha ',\Theta ',\Theta _p'\}$
is the vector of perturbed quantities, c.c. means complex conjugate, and
$\xi ^n(Z)$
is the normal mode amplitude varying in the non-homogeneous direction
$Z$
. Furthermore, according to the temporal stability approach,
$k$
is the real wavenumber, and
$\lambda =\lambda _r+i\lambda _i$
, where
$\lambda _r$
is the temporal growth rate of the perturbation and
$\lambda _i$
is its oscillation frequency. Doing so, leads to their governing equations, namely for the cartesian components
$x$
and
$z$
(subscripts),
\begin{align}&\qquad \lambda \alpha ^n = -\left (\alpha _0 DW^n_{z} + W^n_{z}{D\alpha _0}\right ) - \alpha _0 \dot {\iota } k W^n_{{x}} - \left (W_0 D\alpha ^n + \alpha ^n {D W_0}\right ), \end{align}
\begin{align}& \lambda (D^2 -k^2) {U^n_z} = Pr (D^2 - k^2)^2 U^n_z - Pr Ra k^2 \Theta ^n + \frac {\alpha _0}{2} (\beta -1) \lambda (D^2 -k^2) U_z^n \notag\\ &\qquad\qquad\qquad\quad + {\frac {D\alpha _0 }{2} (\beta -1)\lambda D U_z^n - D\alpha _0{6 Pr (3 -\beta ) }{\Phi ^{-2}} \left(DU_z^n + \dot \iota k W^n_x \right)}\notag \\ &\qquad\qquad\qquad - \alpha _0 {6 Pr (3 -\beta ) }{\Phi ^{-2}}\big(\dot \iota k D W_x^n + (D^2-k^2) U_z^n +k^2 W_z^n\big)\notag\\ &\qquad\qquad\qquad - \big[\tfrac {1}{2}(\beta -1){Ga Pr^{2}} + 6 Pr (3 -\beta ) {W_0}{\Phi ^{-2}} \big]k^2 \alpha ^n, \end{align}
\begin{align}&\qquad \lambda W_z^n + W_0 {D W_z^n} + W_z^n {DW_0} = \beta \lambda {U^n_z} + 12 Pr \beta \Phi ^{-2}{\left(U_z^n-W_z^n\right)}, \end{align}
\begin{align}&\qquad\qquad \lambda W_x^n + W_0 ({DW_x^n}) = \beta \lambda {U^n_x} + 12 Pr \beta \Phi ^{-2}{\left(U_x^n-W_x^n\right)}, \end{align}
\begin{align}& {(1-\alpha _0) } \left [\lambda {\Theta ^n} + U^n_z D \Theta _0\right ] = \alpha ^n \lambda \Theta _0 + (D^2 - k^2) \Theta ^n - 12 \alpha _0 \Phi ^{-2}{\big(\Theta ^n - \Theta _{p}^n\big)}\notag\\ &\qquad\qquad\qquad\qquad\qquad\qquad - 12 \alpha ^n \Phi ^{-2}{\big(\Theta _0^n - \Theta _{p0}^n\big)}, \end{align}
\begin{align}&\qquad\qquad \lambda \Theta ^n_p + W_0 D\Theta ^n_p + W_z^n D \Theta _{p0}^n = {12}({E}{\Phi ^2})^{-1}{\big(\Theta ^n - \Theta ^n_p\big)}, \end{align}
obtained by taking the double curl of the fluid momentum equation to eliminate the pressure, and using the incompressibility relation in order to eliminate the
$x\hbox{-}$
component of the velocity. Their boundary conditions, obtained by similar means from (2.14) and (2.15), are given by
\begin{align} {U_z^n} = {DU_z^n} = 0, \quad \Theta ^n = 0 \quad \text{at} \quad Z= \pm {1}/{2}, \end{align}
\begin{align} \alpha ^n=0, \quad {\boldsymbol{W^n}} = 0, \quad \Theta ^n_{p} = 0\quad \text{at} \quad Z = Z^*. \end{align}
Equations (3.6)–(3.11) are solved numerically by employing the shooting method, where the critical conditions (
$Ra_c, \lambda _c, k_c$
) are determined from the sensitivity of
$Ra$
with respect to
$k$
(see Alves et al. (Reference Alves, Hirata, Schuabb and Barletta2019) for details). In order to verify the accuracy of these results, a matrix-forming approach is also employed. It adopts a fourth-order finite-difference discretisation to transform the differential system of equations into a generalised algebraic eigenvalue problem. The latter is solved numerically by the Arnoldi method with a shift-and-invert spectral transformation (see Souza, Freitas & Alves (Reference de Souza, Freitas and de B. Alves2021) for details). Although not shown here, the excellent agreement between both approaches verifies the accuracy of the numerical results presented in the next section.
4. Results
Since the present study focuses on the influence of particle inlet velocity on the model stability, the following parameters are kept constant throughout the analysis:
$Ga = 12\times 10^{12}$
,
$Pr=5$
,
$E = 5\times 10^{-3}$
,
$\Phi = 10^{-2}$
. Heavy particles are injected from the top with the cold wall temperature
$\Theta _p^* = 0$
, while light particles are injected from the bottom with the hot wall temperature
$\Theta _p^* = 1$
. Furthermore, the inlet flux was varied with respect to a reference value
$\mathcal{J} = \mathcal{J}_0 = 533.3$
. Finally, the thermal feedback of the particles to the fluid is assumed negligible in the present work since
$E\ll 1$
.
4.1. Time-asymptotic analysis
Variation with the density ratio
$\beta$
of the critical Rayleigh number
$Ra_c$
(a) and wave vector
$k_c$
(b). Calculations for various inlet velocities
$W^*$
are shown, from sub-terminal
$\lt W_T$
to super-terminal
$\gt W_T$
. For comparison, the values
$Ra_c,k_c$
corresponding to the single-phase Rayleigh–Bénard (RB) system are drawn as horizontal lines. The particle inlet flux is
$\mathcal{J} = \mathcal{J}_0$
.
The first thing to notice is that the linear onset of instability is characterised by a pitchfork bifurcation in the parameter space explored here. In other words, the oscillatory frequency
$\lambda _i$
is always zero. Second, the effects of the modified density ratio on the critical thresholds under different particle inlet velocity are presented in figure 2. For heavy particles sedimenting from the top (
$\beta \lt 1$
), an increase in the inlet velocity leads to a decrease in the critical Rayleigh number, thereby destabilising the system. Conversely, for light particles injected from the bottom (
$\beta \gt 1$
), the trend is reversed, with higher inlet velocities promoting stabilisation. The most unstable cases occur for sub-terminal particle velocities. Notably, for
$\beta \gt 2.8$
and
$W^*/W_T=0.5$
, the critical Rayleigh number for the pRB model falls below that of the single-phase RB system, indicating that a weakly injected particulate phase can promote instability in the case of light particles and bubbles. It should be noted that the case of neutrally buoyant particles (
$\beta =1$
) is singular, since one would have to inject an infinite amount of particles to keep the volumetric particulate flux constant. As illustrated in figure 2(b), the critical wavenumber slightly decreases with increasing inlet particle velocity for both heavy and light particles, though this decrease is more pronounced in the latter case.
Figure 3(a–d) shows the effect of the particle inlet velocity on the critical thresholds under different inlet volumetric flux intensities for selected values of the modified density ratio. For the case of bubbles (
$\beta =3$
), increasing the inlet flux is destabilising for low inlet velocities but stabilising otherwise, since all three curves intersect at
$W^*/W_T\sim 0.6$
. A similar trend is observed for all light particle cases, though the value of
$W^*/W_T$
where the intersection occurs decreases as
$\beta$
decreases. For
$\beta \lesssim 2.5$
, no intersection is observed and the increase in the particle flux always plays a stabilising role. It is also worth pointing out that, for
$\beta =3$
and
$\mathcal{J}=1.5 \mathcal{J}_0$
, the critical Rayleigh number becomes negative for low enough values of
$W^*/W_T$
. This means that bubbles can trigger convective instabilities even when the pRB model is heated from above. As discussed by Nakamura et al. (Reference Nakamura, Yoshikawa, Tasaka and Murai2020), this is due to a potentially unstable density stratification within the liquid–gas mixture near the bottom wall. As the bubble velocity increases to its terminal velocity, the volume fraction in this region decreases from
$\mathcal{J} /W^*$
to
$\mathcal{J} /W_T$
, increasing the possibility of a locally unstable density gradient.
Critical Rayleigh number as a function of the inlet velocity for three different particle fluxes. The horizontal dashed lines correspond to the single-phase RB threshold.
Same as figure 3, but for the critical wavenumber.
For heavy particles, e.g. figures 3(a) and 3(b), increasing the particle inlet flux is always stabilising, but the effect is more pronounced at low inlet velocities. Additional simulations were conducted for heavy particles, exploring higher values of both particle flux and inlet velocity, within the validity range of the present model. In all cases, the critical Rayleigh number remained above the RB threshold, providing strong evidence that it is not possible to destabilise the RB model by adding heavy particles. Furthermore, figure 4(a–d) shows that the wavenumber dependence on both particle flux and inlet velocity is relatively weaker. This means that the size of the convective cells at the onset of convection are not significantly different from their single-phase RB counterparts. Such a difference becomes more pronounced, however, for very slow-rising bubbles at high volumetric inlet rates, as illustrated by the red line in figure 4(d).
Streamlines and particle concentration at the onset of convection are plotted in figure 5, through the eigenvectors of the perturbed fluid velocity
$\textbf{U}'$
and particle distribution
$\alpha '$
, for both heavy and light particles. Heavy particles (figure 5a–c
) accumulate near the bottom. For sub-terminal (super-terminal) inlet velocities, a higher particle volume fraction is observed in the downwelling (upwelling) plumes. The particular case where
$W^*=W_T$
is not shown for conciseness as it is already known that it possess a uniform particle distribution (Raza et al. Reference Raza, Hirata and Calzavarini2024). The opposite behaviour is observed for light particles (figure 5d–f
). They accumulate near the top, with upwelling (downwelling) plumes favouring higher concentrations for sub-terminal (super-terminal) velocities. Finally, when comparing figures 5(c) and 5(f), it can be observed that heavy particles (
$\beta =0.5$
) tend to accumulate closer to the walls than light ones (
$\beta =3$
) at super-terminal velocities (
$W^*/W_T=1.2$
).
Streamlines of fluid velocity field and colourmap of the particle volume fraction at the onset of convection for
$\beta = 0.5$
(a–c) and
$\beta = 3$
(d–f) and increasing inlet velocities from left to right. For a better comparison, the particle volume fraction
$\alpha '$
is normalised with respect to the base volume fraction
$\alpha _0$
for
$\beta =3$
and
$W^*=W_T$
particles. The particle inlet flux is
$\mathcal{J} = \mathcal{J}_0$
.
4.2. Energy budget analysis
An energy budget analysis is performed in order to gain further insights into the physical mechanisms triggering instabilities in the pRB system. In order to derive the evolution equation for the perturbation kinetic energy, we take the product of the momentum equation (3.7) and the complex conjugate of the vertical fluid velocity,
$\bar {U}_z^n$
, where the overbar denotes complex conjugate. By integrating the resulting equation over the entire domain, the following relationship for the spatially averaged perturbation kinetic energy rate is obtained:
\begin{equation} \lambda \ e^K= e^{\Theta } + e^{V} + e^{\alpha _0W} + e^{\alpha W_0}+e^{\alpha Ga}, \end{equation}
Components of the kinetic energy budget at the neutral stability condition versus the particle inlet velocity: thermal buoyancy injection power
$E^{\Theta }$
, viscous dissipation rate
$E^{V}$
, particle’s feedback due to the base particle concentration, base particle velocity and particle buoyancy
$E^{\alpha _0 W}$
,
$E^{\alpha W_0}$
,
$E^{\alpha Ga}$
. Results obtained for different particle inlet fluxes:
$0.5\mathcal{J}_0$
(dotted lines),
$\mathcal{J}_0$
(dashed lines),
$1.5\mathcal{J}_0$
(solid lines).
with each term in this equation being defined as
\begin{align}&\qquad\qquad e^K = \left (1 - \dfrac {\alpha _0}{2} (\beta -1) \right ) {\int _{-1/2}^{1/2}} \mathrm{Re}\big[(D^2 -k^2) (U_z^n) \bar {U}_z^n \big]{\rm d}Z\nonumber \\ &\qquad\qquad\quad{\quad +\dfrac {(\beta -1)}{2}{\int _{-1/2}^{1/2}} D\alpha _0\mathrm{Re}\left [D (U_z^n) \bar {U}_z^n \right ]{\rm d}Z}, \end{align}
\begin{align}&\qquad\qquad\qquad\quad\,\, e^{\Theta } = Pr Ra {\int _{-1/2}^{1/2}} \mathrm{Re}\big[(\Theta ^n k^2) \bar {U}_z^n\big]{\rm d}Z, \end{align}
\begin{align}&\qquad\qquad\qquad e^{V} = -Pr {\int _{-1/2}^{1/2}} \mathrm{Re}\big[(D^2 -k^2)^2 (U_z^n) \bar {U}_z^n \big ]{\rm d}Z, \end{align}
\begin{align} e^{\alpha _0 W} &= \dfrac {6\, Pr (3 - \beta )}{\Phi ^2}\bigg [ {\int _{-1/2}^{1/2}} \alpha _0 \, \mathrm{Re} \big[ (D^2 - k^2) (U_z^n) \bar {U}_z^n + (\dot {\iota } k D W_x^n + k^2 W_z^n) \bar {U}_z^n \big]{\rm d}Z \nonumber \\ &{\quad +{\int _{-1/2}^{1/2}} D\alpha _0 \mathrm{Re} \big[ (\dot {\iota } k W_x^n)\bar {U}_z^n + D (U_z^n) \bar {U}_z^n \big]} \ {\rm d}Z\bigg ], \end{align}
\begin{align}&\qquad\qquad\quad e^{\alpha W_0} = \dfrac {6 Pr (3 -\beta )}{\Phi ^2} {\int _{-1/2}^{1/2}} W_0 \mathrm{Re}\left [ (\alpha ^n k^2) \bar {U}_z^n \right ]{\rm d}Z, \end{align}
\begin{align}&\qquad\qquad\quad e^{\alpha Ga} = {\dfrac {(\beta -1) {Ga Pr^2}}{2}}{\int _{-1/2}^{1/2}} \mathrm{Re}\left [ (\alpha ^n k^2) \bar {U}_z^n \right ]{\rm d}Z,\end{align}
where
$e^{\Theta }$
is the power of the thermal buoyancy force, and
$e^{V}$
indicates the rate of viscous energy dissipation. Furthermore,
$e^{\alpha _0 W}$
,
$e^{\alpha W_0}$
and
$e^{\alpha Ga}$
collectively represent the particle feedback, i.e. the parts due to the base particle concentration, base particle velocity and particle buoyancy, respectively. The first two terms represent the drag force, while the latter represents the Archimedes force. The superscript
$\alpha$
denotes the terms arising from the non-homogeneous distribution of particles, which are absent if
$W^*=W_T$
as in Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021) and Raza et al. (Reference Raza, Hirata and Calzavarini2024). By normalising these contributions using the absolute value of the dissipation rate, we obtain at neutral condition
$(\lambda = 0)$
\begin{equation} E^{\Theta } + E^{\alpha _0W} + E^{\alpha W_0}+ E^{\alpha Ga} = 1. \end{equation}
Figure 6 shows the aforementioned normalised kinetic energy rates, collectively denoted as
$E_K$
, versus the particle inlet velocity for different particle fluxes
$ \mathcal{J}$
and density ratios
$ \beta$
. Positive (negative) contributions indicate destabilisation (stabilisation) of the steady state. Heavy-particle cases (figure 6
a,b) show a competition between the destabilising thermal buoyancy and the stabilising particle drag feedback. Hence, heavy-particle action is purely dissipative as their friction opposes the thermal buoyancy. However, the energy injection due to particle buoyancy
$ E^{\alpha Ga}$
is always negligible in the present study. This is at odds with the light-particle cases (figure 6
c,d) where the particle buoyancy term turns out to be important. Here
$ E^{\alpha Ga}$
is destabilising when the injection velocity is smaller than the terminal velocity, while the opposite is true otherwise. This highlights the importance of the direction of particle motion, namely ascending or descending, as well as the sign of their acceleration, which determines whether the particle concentration near the injection wall is diluted or intensified. We note that in the special case of bubbles (figure 6
d), where the drag terms
$ E^{\alpha _0 W}$
and
$ E^{\alpha W_0}$
vanish by definition according to (4.5)–(4.6), the competition is solely between the thermal and particle buoyancies.
5. Concluding remarks
The present investigation highlighted how the particle inlet velocity influences the linear stability of the pRB system. Increasing the inlet velocity while maintaining a constant particle flux destabilises the system for heavy particles but progressively stabilises it for light particles. While the general features of particle accumulation persist across the range of inlet velocities, the spatial localisation of accumulation shifts between upwelling and downwelling regions as the injection speed transitions from sub-terminal to super-terminal values. When the injection velocity matches the terminal velocity, the linearised system dynamics does not support the accumulation of particles. In that case, particle accumulation could arise only through nonlinear interactions as shown by Srinivas & Tomar (Reference Srinivas and Tomar2025). In this study, the thermal coupling was deliberately minimised in order to isolate the effects of the mechanical coupling introduced by the injection velocity. This was achieved by maintaining
$E \ll 1$
and setting the particle inlet temperature equal to the inlet-wall temperature. Nevertheless, particle thermal inertia and injection temperature offer alternative mechanisms for influencing the flow onset within this model. Preliminary investigations in this direction (Raza, Calzavarini & Hirata Reference Raza, Calzavarini and Hirata2025) indicate that increasing
$E$
consistently enhances system stability when particles are heavier than the fluid, regardless of the volumetric particle flux, injection velocity or injection temperature. In contrast, for particles lighter than the fluid, the effects are more complex: the influence of
$E$
can be either stabilising or destabilising, depending on the inlet velocity and particulate volumetric flux. This suggest that the effects of momentum coupling dominate over the thermal ones. However, we shall note that the parameter
$E$
depends on the relative density between the particles and the fluid, in fact
$E=(c_{Pp}/c_P)(3-\beta )/(2\beta )$
, implying that when the particles are very light, the thermal coupling decreases and eventually becomes negligible in the bubble limit. The physical mechanisms underlying the behaviour of the thermal coupling remain to be fully understood and will be the focus of future work.
Finally, an essential direction for future investigation of this problem is the impact of particle size. This can be addressed only in a very limited way within the framework of an Eulerian model, such as that considered here – for example, by including Faxén corrections in the particle equations. However, it can be more effectively approached through numerical studies with resolved particles, such as those recently proposed by Chen & Prosperetti (Reference Chen and Prosperetti2024).
Declaration of interests
The authors report no conflict of interest.
