2. Problem statement and mathematical model

The physical model and domain with schematic representation, as shown in figure 1, consists of a continuous fluid domain $\varOmega =\left \{(x^*,z^*)\in \mathbb{R}\times (-H/2, H/2)\right \}$ , along with a dispersed particulate phase. Here, the bottom and the top surfaces are treated as isothermal walls maintained at temperatures $T_h$ and $T_c$ , respectively. The present study considers the monodispersion of very tiny spherical particles at relatively small volume fractions ( $\leqslant 10^{-3}$ ). There are two alternative methods, called Euler–Euler and Euler–Lagrange formulations, available for the particles in the millimetre range for which the homogeneous equilibrium model fails. However, the choice of a particular method depends on various length scales, such as particle size $d_{\!p}$ , the smallest flow length scale $\varDelta$ , inter-particle separation $\lambda$ , and particle number density $n$ . When $d_{\!p}\ll \varDelta \sim \lambda$ and $n\varDelta ^3\gg 1$ , the Euler–Lagrange formulation is preferred, whereas if $d_{\!p}\ll \lambda \ll \varDelta$ and $n\Delta ^3\gg 1$ , then the Euler–Euler, also called the two-fluid formulation, is preferred. The inter-particle separation $\lambda$ is related to the non-dimensional particle size $\delta =(d_{\!p}/H)$ and volume fraction as $(\lambda /H)\sim (\delta /\varPhi _0^{1/3})$ , and, the condition becomes $\delta \ll (\delta /\varPhi _0^{1/3})\ll (\varDelta /H)$ . Here, in the present work, we consider the particle volume fraction $\sim 10^{-3}$ , thus the first inequality, $\delta \ll (\delta /\varPhi _0^{1/3})$ , is satisfied. However, the second inequality is satisfied only when the particle size is $\delta \ll (\varDelta /H)/\varPhi _0^{1/3}$ . Therefore, when $(\varDelta /H)\sim 10^{-1}$ and $\varPhi _0\sim 10^{-3}$ , the particle size is $\delta \ll 1$ . Hence in the present work, we choose $\delta =0.01,0.05$ . We adopt the two-fluid model given by Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021), in which the particles are introduced steadily and uniformly at the top surface at their terminal velocity with a fixed temperature. The momentum equation for the particulate phase is given by the volume-averaged Maxey–Riley–Gatignol equation (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983), and the thermal energy equation ensures the balance between the rate of change of the sensible heat of the particles and the rate of convective heat transfer between particulate and fluid. The dimensional independent variables $ (x^*,z^*,t^* )$ and dependent variables $ (u_x^*,u_z^*,p^*,T^*,v_x^*,v_z^*,T_{\!p}^* )$ are non-dimensionalised as follows:

(2.1) \begin{equation} \left .\begin{array}{l} \left (x,z,t\right )=\left (x^*/H,z^*/H,t^*U/H\right )\!,\\[8pt]\left (u_x,u_z,p,\theta \right )=\left (u_x^*/U,u_z^*/U,p^*/\left(\rho _{\!f}U^2\right)\!,(T-T_c)/\Delta T\right )\!,\\[8pt](v_x,v_z,\theta _{\!p})=\left (v_x^*/U,v_z^*/U,(T_{\!p}-T_c)/\Delta T\right )\!, \end{array}\!\!\right \} \end{equation}

where $x$ , $z$ and $t$ are the non-dimensional horizontal coordinate, vertical coordinate and time, respectively. Furthermore, $u_x$ , $u_z$ , $p$ and $\theta$ are the fluid horizontal velocity component, vertical velocity component, pressure and temperature, respectively. Similarly, $v_x$ , $v_z$ and $\theta _{\!p}$ are dimensionless particulate phase horizontal velocity component, vertical velocity component and temperature, respectively. Here, we use the distance between the two horizontal surfaces $H$ as a length scale, fluid free-fall velocity $U=\sqrt {g\beta\, \Delta \textit{TH}}$ as a velocity scale, and the temperature difference $\Delta T=T_h-T_c$ as a temperature scale, where $\rho _{\!f}$ , $g$ and $\beta$ are fluid density, the acceleration due to gravity and the volume expansion coefficient of fluid, respectively. The resulting governing equations in non-dimensional form are given as

Figure 1. Schematic of the particle-laden Rayleigh–Bénard convection.

(2.2) \begin{align} \frac {\partial u_x}{\partial x}+\frac {\partial u_z}{\partial z}&=0, \end{align}

(2.3) \begin{align} \frac {\partial u_x}{\partial t}+u_x\frac {\partial u_x}{\partial x}+u_z\frac {\partial u_x}{\partial z}&=-\frac {\partial p}{\partial x} + \sqrt {\frac {\textit{Pr}}{\textit{Ra}}}{\nabla}^{2}u_x-\frac {{R}\phi }{{\textit{St}}_m}{\left (u_x-v_x\right )}, \end{align}

(2.4) \begin{align} \frac {\partial u_z}{\partial t}+u_x\frac {\partial u_z}{\partial x}+u_z\frac {\partial u_z}{\partial z}&=-\frac {\partial p}{\partial z} + \sqrt {\frac {\textit{Pr}}{\textit{Ra}}}{\nabla}^{2}u_z+\theta -\frac {{R}\phi }{{\textit{St}}_m}{\left (u_z-v_z\right )}, \end{align}

(2.5) \begin{align} \frac {\partial \theta }{\partial t}+u_x\frac {\partial \theta }{\partial x}+u_z\frac {\partial \theta }{\partial z}&=\frac {1}{\sqrt {{{\textit{Ra}}{\textit{Pr}}}}}{\nabla}^2 \theta -\frac {{E}\phi }{{\textit{St}}_{\textit{th}}}{(\theta -\theta _{\!p})}, \end{align}

(2.6) \begin{align} \frac {\partial \phi }{\partial t}+\frac {\partial {\left (\phi v_x\right )}}{\partial x}+\frac {\partial \left (\phi v_z\right )}{\partial z}&=0, \end{align}

(2.7) \begin{align} \frac {\partial v_x}{\partial t}+v_x\frac {\partial v_x}{\partial x}+v_z\frac {\partial v_x}{\partial z}&=\frac {u_x-v_x}{{\textit{St}}_m}, \end{align}

(2.8) \begin{align} \frac {\partial v_z}{\partial t}+v_x\frac {\partial v_z}{\partial x}+v_z\frac {\partial v_z}{\partial z}&=\frac {u_z-v_z}{{\textit{St}}_m}-\frac {v_0}{{\textit{St}}_m}, \end{align}

(2.9) \begin{align} \frac {\partial \theta _{\!p}}{\partial t}+v_x\frac {\partial \theta _{\!p}}{\partial x}+v_z\frac {\partial \theta _{\!p}}{\partial z}&=\frac {\theta -\theta _{\!p}}{{\textit{St}}_{\textit{th}}} \end{align}

for all $ (x,z,t )\in \left \{\mathbb{R}\times (-1/2, 1/2 )\times (0, \infty )\right \}$ . The boundary conditions are

(2.10) \begin{equation} \left .\begin{array}{l} \mbox{at }z=1/2,\quad u_x=u_z=\theta =v_x=0,\mbox{ }v_z=-v_0,\mbox{ }\theta _{\!p}=\varTheta _{\textit{pt}}\mbox{ and }\phi =\varPhi _0,\\[3pt] \mbox{at }z=-1/2,\quad u_x=u_z=0\mbox{ and }\theta =1, \end{array}\!\!\right \} \end{equation}

for all ${ (x,t )}\in \left \{\mathbb{R}\times { (0,\infty )}\right \}$ , where $\phi$ is the particle volume fraction field, and the non-dimensional parameters present in the problem, Rayleigh number ${\textit{Ra}}$ , Prandtl number $ \textit{Pr} $ , density ratio $R$ , specific heat capacity ratio $E$ , mechanical Stokes number $St_m$ , thermal Stokes number $St_{\textit{th}}$ , initial particle temperature $\varTheta _{\textit{pt}}$ , and non-dimensional particle terminal speed $v_0$ , are defined as

(2.11a–h) \begin{align} &\textit{Ra}=g\beta\, \Delta {\textit{TH}}^3/(\alpha _{\!f}\nu _{\!f}), \quad \textit{Pr}=\nu _{\!f}/\alpha _{\!f}, \quad R=\rho _{\!p}/\rho _{\!f}, \nonumber \\&{\textit{St}}_m=\frac {\tau _{\!p}}{\tau _{\!f}}=\frac {{R}\delta ^2}{18{C}_d} \sqrt {\frac {\textit{Ra}}{\textit{Pr}}},\quad {\textit{St}}_{\textit{th}}=\frac {\tau _{\textit{th}}}{\tau _{\!f}}=\frac {{E}\delta ^2}{6{\textit{Nu}}_{\!p}}\sqrt {{{\textit{Ra}}{\textit{Pr}}}}, \nonumber \\&E=R(C_{\textit{pp}}/C_{\textit{pf}}), \quad v_0=\frac {{\textit{Re}}_{\!p}}{\delta }\sqrt {\frac {\textit{Pr}}{\textit{Ra}}}, \quad \varTheta _{\textit{pt}}=\frac {T_{\textit{pt}}-T_c}{T_h-T_c}, \end{align}

where $\alpha _{\!f}=\kappa _{\!f}/\rho _{\!f}C_{\textit{pf}}$ is the fluid thermal diffusion coefficient, $\nu _{\!f}=\mu _{\!f}/\rho _{\!f}$ is the fluid kinematic viscosity, $\kappa _{\!f}$ is the fluid thermal conductivity, $C_{\textit{pf}}$ is the fluid specific heat constant, $\mu _{\!f}$ is fluid dynamic viscosity, $\rho _{\!p}$ is the particle density, and $C_{\textit{pp}}$ is the particle specific heat. Particle diameter $d_{\!p}$ in the non-dimensional form is represented as $\delta =d_{\!p}/H$ . Introduction of particles into the flow leads to three distinct time scales: $\tau _{\!f}=H/U$ is the flow time scale, $\tau _{\!p}=\rho _{\!p}d_{\!p}^2/(18\mu _{\!f}C_d({\textit{Re}}_{\!p}))$ is the particle mechanical relaxation time scale, and $\tau _{\textit{th}}=\rho _{\!p}C_{\textit{pp}}d_{\!p}^2/(6\kappa _{\!f}\,{\textit{Nu}}_{\!p}({\textit{Re}}_{\!p},\textit{Pr}))$ is the particle thermal relaxation time scale. Here, ${C}_d=1+0.15\,{\textit{Re}}_{\!p}^{0.687}$ is the Schiller–Naumann correction factor to the Stokes drag when the particle Reynolds number ${\textit{Re}}_{\!p}=\rho _{\!f}v_\tau d_{\!p}/\mu _{\!f}$ is greater than unity and less than 800 (Clift, Grace & Weber Reference Clift, Grace and Weber2005), $v_\tau ={ (1-\rho _{\!f}/\rho _{\!p} )}\tau _{\!p}g$ is the settling speed of a particle in the quiescent fluid, and ${\textit{Nu}}_{\!p}=2.0+0.6\,{\textit{Re}}_{\!p}^{1/2}\,\textit{Pr}^{1/3}$ is the particle Nusselt number. The last terms on the right-hand sides of (2.3) and (2.4) represent the mechanical two-way coupling between the fluid and particles, whereas the thermal two-way coupling is captured by the right-hand-side last term of (2.5).

3. Linear stability analysis

We assume that the fluid–particle system is at a steady state, with particles settling uniformly with their terminal velocity in a quiescent fluid. Under these conditions, the above governing equations (2.2)–(2.9) reduced to a system of ordinary differential equations, which are represented in operator form as

(3.1) \begin{align} \mathcal{L}_z{\left (P_0,\varTheta _0,v_0,\varPhi _0,{R},{\textit{St}}_m\right )}&=\frac {\textrm{d} P_0}{\textrm{d} z}-\varTheta _0+\frac {{R}\varPhi _0}{{\textit{St}}_m}v_0=0, \\[-12pt] \nonumber \end{align}

(3.2) \begin{align} \mathcal{L}_\theta {\left (\varTheta _0,\varTheta _{p0},\varPhi _0,{\textit{Ra}},{\textit{Pr}},{E}\right )}&=\frac {1}{\sqrt {{{\textit{Ra}}{\textit{Pr}}}}}\frac {\textrm{d}^2\varTheta _0}{\textrm{d}z^2}-\frac {{E}\varPhi _0}{{\textit{St}}_{\textit{th}}}{\left (\varTheta _0-\varTheta _{p0}\right )}=0, \\[-12pt] \nonumber \end{align}

(3.3) \begin{align} \mathcal{L}_{p\theta }{\left (\varTheta _0,\varTheta _{p0},v_0,{\textit{St}}_{\textit{th}}\right )}&=\frac {\varTheta _0-\varTheta _{p0}}{{\textit{St}}_{\textit{th}}}+v_0\frac {\textrm{d} \varTheta _{p0}}{\textrm{d} z}=0 \\[9pt] \nonumber \end{align}

for $z\in { (-1/2,1/2 )}$ with boundary conditions

(3.4) \begin{equation} \left .\begin{array}{l} \mbox{at }z=1/2,\quad\varTheta _0=0\mbox{ and } \varTheta _{p0}=\varTheta _{\textit{pt}},\\[3pt] \mbox{at }z=-1/2,\quad\varTheta _0=1, \end{array}\!\!\right \} \end{equation}

where $P_0$ , $\varTheta _0$ and $\varTheta _{p0}$ are the basic state fluid pressure, basic state fluid temperature and basic state particle temperature field, respectively, and $\varPhi _0$ is the initial uniform particle volume fraction. Using the above boundary conditions, the solution to (3.1)–(3.3) yields the base state given by

(3.5) \begin{align} \varTheta _0&=\varTheta _{\textit{pt}}+\mathscr{A}\big [1-{\left (1-{\textit{lm}}_1\right )}\,e^{-m_1(1/2-z)}\big ]+\mathscr{B}\big [1-{\left (1-{\textit{lm}}_2\right )}e^{-m_2(1/2-z)}\big ], \\[-9pt] \nonumber \end{align}

(3.6) \begin{align} \varTheta _{p0}&=\varTheta _{\textit{pt}}+\mathscr{A}\big [1-{\textrm{e}}^{-m_1(1/2-z)}\big ]+\mathscr{B}\big [1-{\textrm{e}}^{-m_2(1/2-z)}\big ], \\[-9pt] \nonumber \end{align}

(3.7) \begin{align} \frac {\textrm{d} P_0}{\textrm{d} z}&=\varTheta _{\textit{pt}}+\mathscr{A}\big [1-{\left (1-{\textit{lm}}_1\right )}\,e^{-m_1(1/2-z)}\big ] + \mathscr{B}\big [1-{\left (1-{\textit{lm}}_2\right )}e^{-m_2(1/2-z)}\big ]- \frac {{R}\varPhi _0}{{\textit{St}}_m}v_0, \\[9pt] \nonumber \end{align}

where $\mathscr{A}$ and $\mathscr{B}$ are integration constants given by

(3.8) \begin{align} \mathscr{A}&=-\frac {(1-\varTheta _{\textit{pt}}){\textit{lm}}_2+\varTheta _{\textit{pt}}\left [1-{\left (1-{\textit{lm}}_2\right )}\,e^{-m_2}\right ]}{{\textit{lm}}_1\left [1-{\left (1-{\textit{lm}}_2\right )}\,e^{-m_2}\right ]-{\textit{lm}}_2\left [1-{\left (1-{\textit{lm}}_1\right )}\,e^{-m_1}\right ]}, \end{align}

(3.9) \begin{align} \mathscr{B}&=+\frac {(1-\varTheta _{\textit{pt}}){\textit{lm}}_1+\varTheta _{\textit{pt}}\left [1-{\left (1-{\textit{lm}}_1\right )}\,e^{-m_1}\right ]}{{\textit{lm}}_1\left [1-{\left (1-{\textit{lm}}_2\right )}\,e^{-m_2}\right ]-{\textit{lm}}_2\left [1-{\left (1-{\textit{lm}}_1\right )}\,e^{-m_1}\right ]}. \end{align}

Here, $m_1$ and $m_2$ are given by

(3.10a,b) \begin{align} m_1=\frac {1}{2l}\left [1+\sqrt {1+4{\left (\frac {l}{{L}}\right )}^2}\right ]\!,\quad m_2=\frac {1}{2l}\left [1-\sqrt {1+4{\left (\frac {l}{{L}}\right )}^2}\right ]\!, \end{align}

where $L$ is the spatial non-dimensional length scale over which the effect of a particle on the surrounding fluid temperature is significant (Prakhar & Prosperetti Reference Prakhar and Prosperetti2021), given by

(3.11) \begin{equation} {L}={\left (\frac {{\textit{St}}_{\textit{th}}}{{E}\varPhi _0\sqrt {{{\textit{Ra}}{\textit{Pr}}}}}\right )}^{1/2}=\frac {\delta }{\sqrt {6\varPhi _0\,{\textit{Nu}}_{\!p}}}. \end{equation}

Here, $l=v_0{\textit{St}}_{\textit{th}}$ is the non-dimensional length scale that the particle must traverse for its temperature to be locally equal to that of the fluid (Prakhar & Prosperetti Reference Prakhar and Prosperetti2021), related to other parameters as

(3.12) \begin{equation} l=\frac {{\textit{St}}_{\textit{th}}\,{\textit{Re}}_{\!p}}{\delta }\sqrt {\frac {\textit{Pr}}{\textit{Ra}}}=\frac {{E}\delta\, {\textit{Pr}}\,{\textit{Re}}_{\!p}}{6\,{\textit{Nu}}_{\!p}}. \end{equation}

The classical normal mode analysis (Drazin & Reid Reference Drazin and Reid2004) is performed to examine the stability of the basic flow mentioned above. In general, the linear stability analysis is carried out by decomposing all the dependent variables into a steady basic state and the respective infinitesimal perturbations (represented by a superscript prime):

(3.13) \begin{align} {\big (u_x,u_z,\theta ,v_x,v_z,\theta _{\!p},p,\phi \big )}={}&{\big (0,0,\varTheta _0(z),0,-v_0,\varTheta _{p0}(z),P_{0}(z),\varPhi _0\big )}\nonumber\\&{}+{\big (u_x^{\prime},u_z^{\prime},\theta ',v_x^{\prime},v_z^{\prime},\theta _{\!p}^{\prime},p',\phi '\big )}. \end{align}

After neglecting the higher-order terms and retaining only the first-order terms in $ \{u_x^{\prime},u_z^{\prime},\theta ',v_x^{\prime},v_z^{\prime},\theta _{\!p}^{\prime},p',\phi ' \}$ , we obtain

(3.14) \begin{align} \frac {\partial u_x^{\prime}}{\partial x}+\frac {\partial u_z^{\prime}}{\partial z}&=0, \end{align}

(3.15) \begin{align} \frac {\partial u_x^{\prime}}{\partial t}&=-\frac {\partial p'}{\partial x}+\sqrt {\frac {\textit{Pr}}{\textit{Ra}}}{\nabla}^2{u_x^{\prime}}-\frac {{R}\varPhi _0}{{\textit{St}}_m}{\left (u_x^{\prime}-v_x^{\prime}\right )}, \end{align}

(3.16) \begin{align} \frac {\partial u_z^{\prime}}{\partial t}&=-\frac {\partial p'}{\partial z}+\sqrt {\frac {\textit{Pr}}{\textit{Ra}}}{\nabla}^2{u_z^{\prime}}+\theta '-\frac {{R}\varPhi _0}{{\textit{St}}_m}{\left (u_z^{\prime}-v_z^{\prime}\right )}-\frac {{R}v_0}{{\textit{St}}_m}\phi ', \end{align}

(3.17) \begin{align} \frac {\partial \theta '}{\partial t}&=\frac {1}{\sqrt {{{\textit{Ra}}{\textit{Pr}}}}}{\nabla}^2\theta '-u_z^{\prime}\frac {\textrm{d} \varTheta _0}{\textrm{d} z}-\frac {{E}\varPhi _0}{{\textit{St}}_{\textit{th}}}{\big (\theta '-\theta _{\!p}^{\prime}\big )}-\frac {{E}{\big (\varTheta _0-\varTheta _{p0}\big )}}{{\textit{St}}_{\textit{th}}}\phi ', \end{align}

(3.18) \begin{align} \frac {\partial \phi '}{\partial t}&=v_0\frac {\partial \phi '}{\partial z}-\varPhi _0{\left (\frac {\partial v_x^{\prime}}{\partial x}+\frac {\partial v_z^{\prime}}{\partial z}\right )}, \end{align}

(3.19) \begin{align} \frac {\partial v_x^{\prime}}{\partial t}&=v_0\frac {\partial v_x^{\prime}}{\partial z}+\frac {u_x^{\prime}-v_x^{\prime}}{{\textit{St}}_m}, \end{align}

(3.20) \begin{align} \frac {\partial v_z^{\prime}}{\partial t}&=v_0\frac {\partial v_z^{\prime}}{\partial z}+\frac {u_z^{\prime}-v_z^{\prime}}{{\textit{St}}_m}, \end{align}

(3.21) \begin{align} \frac {\partial \theta _{\!p}^{\prime}}{\partial t}&=v_0\frac {\partial \theta _{\!p}^{\prime}}{\partial z}+\frac {\theta '-\theta _{\!p}^{\prime}}{{\textit{St}}_{\textit{th}}}. \\[12pt] \nonumber \end{align}

The above system of equations is linear homogeneous and has the coefficient functions of spatial variables only, not dependent on time. Hence they satisfy the solution in the normal mode form given by

(3.22) \begin{align} &{\big (u_x^{\prime},u_z^{\prime},\theta ',v_x^{\prime},v_z^{\prime},\theta _{\!p}^{\prime},p',\phi '\big )}^{\textrm T}\nonumber \\&\quad=e^{{\textrm i}k{ (x-ct )}}{\big (\hat {u}_x(z),\hat {u}_z(z),\hat {\theta }(z),\hat {v}_x(z),\hat {v}_z(z),\hat {\theta }_{\!p}(z),\hat {p}(z),\hat {\phi }(z)\big )}^{\textrm T}, \end{align}

where $k$ is the wavenumber, and $c=c_r+\textrm{i}c_r$ is the complex wave speed corresponding to $k$ . The sign of $c_i$ determines the growth or decay of the disturbance. That is, as $c_i\lt 0$ or $c_i=0$ or $c_i\gt 0$ , the flow is stable or neutrally stable or unstable, respectively. Substitution of (3.22) in (3.14)–(3.21) leads to linearised disturbance equations in operator form given in the § A.1. The boundary conditions for corresponding equations are given by

(3.23) \begin{equation} \left .\begin{array}{l} \mbox{at}\ z=1/2,\quad\hat {u}_x=\hat {u}_z=\hat {\theta }=\hat {v}_x=\hat {v}_z=\hat {\theta }_{\!p}=\hat {\phi }=0,\\ \mbox{at}\ z=-1/2,\quad\hat {u}_x=\hat {u}_z=\hat {\theta }=0. \end{array}\right \} \end{equation}

The system of linear equations (A1)–(A8) together with the boundary conditions gives rise to a generalised eigenvalue problem given by

(3.24) \begin{equation} \mathbb{A}\boldsymbol{q}=c\mathbb{B}\boldsymbol{q}, \end{equation}

where $\boldsymbol{q}={ (\hat {u}_x,\hat {u}_z,\hat {\theta },\hat {v}_x,\hat {v}_z,\hat {\theta }_{\!p},\hat {p},\hat {\phi } )}$ is the eigenfunction corresponding to the eigenvalue $c$ , and $ \mathbb{A},\mathbb{B} $ are the square complex matrices. The fundamental disturbance is given by $\boldsymbol{q}\,e^{{\textrm i}k{ (x-c_rt )}}$ , where $\boldsymbol{q}$ is the eigenfunction related to the least stable eigenvalue, and $c_r$ is the corresponding wave speed. (Several parts of the weakly nonlinear analysis will employ the frequent use of the fundamental disturbance.)

The bifurcation point, also known as the critical point $(k, Ra)$ , and the shape of the emerging disturbances are both determined using the linear stability theory. It offers no information regarding the actual magnitude of the disturbances (amplitude) away from the critical value. A weakly nonlinear stability study is necessary to examine the amplitude of such disturbances. The outcomes of the linear stability analysis are also necessary for a nonlinear analysis.

There have been two studies (Prakhar & Prosperetti Reference Prakhar and Prosperetti2021; Raza et al. Reference Raza, Hirata and Calzavarini2024) on the linear stability analysis of the present problem in the literature. Using a two-fluid model, Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021) studied the stability threshold of particle-laden Rayleigh–Bénard convection when the particles are much denser than the underlying fluid. Raza et al. (Reference Raza, Hirata and Calzavarini2024) extended the two-fluid model by Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021) to lighter particles such as bubbles by adding an added-mass term. Using dimensional analysis, we can show that eight independent dimensionless parameters exist in this problem, such as Rayleigh number ${\textit{Ra}}$ , Prandtl number $ \textit{Pr} $ , density ratio $R$ , heat capacity ratio $E$ , particle Reynolds number ${\textit{Re}}_{\!p}$ , particle diameter $\delta$ , particle injection temperature $\varTheta _{\textit{pt}}$ , and initial particle volume fraction $\varPhi _0$ . We fix the dimensionless parameters such as Prandtl number $ \textit{Pr}=0.71$ , density ratio $R=800$ , and specific heat capacity ratio $E=3385$ , such that the fluid–particle system represents the water droplets suspended in dry air.

We present the linear stability results for the air–water droplet system. The effect of initial particle volume fraction on the critical Rayleigh number ${\textit{Ra}}_c$ and the critical wavenumber $k_c$ for different particle Reynolds numbers ${\textit{Re}}_{\!p}$ and the particle sizes $\delta$ is shown in figure 2(a,b). The critical Rayleigh number increases with $\varPhi _0$ for all ${\textit{Re}}_{\!p}$ and $\delta$ . This stabilisation of the flow under the presence of particles is due to the dissipative nature of the mechanical and thermal two-way coupling source terms in (2.3)–(2.4) and (2.5), respectively. The strength of the mechanical source terms is proportional to $({R}\varPhi _0/{\textit{St}}_m)\sim \varPhi _0{C}_d/\delta ^2$ , independent of the density ratio $R$ . Similarly, the strength of the thermal source term is proportional to ${E}\varPhi _0/{\textit{St}}_{\textit{th}}\sim \varPhi _0\,{\textit{Nu}}_{\!p}/\delta ^2$ , which is independent of the heat capacity ratio $E$ . Hence with an increase in the $\varPhi _0$ value (decrease in $\delta$ value) in either source term, the dissipation increases and thereby increases the stability (increase in ${\textit{Ra}}_c$ value). However, the increase in ${\textit{Ra}}_c$ with $\varPhi _0$ is more significant for smaller particles ( $\delta =0.01$ ), due to the strong $1/\delta ^2$ dependence of the two-way coupling source terms, than for the large particles ( $\delta =0.05$ ). Moreover, the weak dependency of the mechanical and thermal source terms on particle Reynolds number ${\textit{Re}}_{\!p}$ through drag force coefficient ${C}_d$ and the Nusselt number ${\textit{Nu}}_{\!p}$ explains the small increase in ${\textit{Ra}}_c$ with the increase in ${\textit{Re}}_{\!p}$ . Similarly, from figure 2(b), it is clear that the critical wavenumber $k_c$ increases with an increase in $\varPhi _0$ for all ${\textit{Re}}_{\!p}$ and $\delta$ . The plausible explanation for this can be obtained using the inter-particle distance $\lambda /H\sim (1-\varPhi _0)^{1/3}\delta /\varPhi _0^{1/3}\approx \delta /\varPhi _0^{1/3}$ for the dilute suspensions $(\varPhi _0\ll 1)$ (Prakhar & Prosperetti Reference Prakhar and Prosperetti2021). Hence the critical wavenumber can be expected to vary as $k_c\sim {\varPhi _0^{1/3}}/{\delta }$ , which explains the increase of $k_c$ with an increase in $\varPhi _0$ for all ${\textit{Re}}_{\!p}$ and $\delta$ values. We note that the increase in $k_c$ is significant for the small particles ( $\delta =0.01$ ) with small Reynolds number ${\textit{Re}}_{\!p}=1$ rather than for the large particles ( $\delta =0.05$ ) with large Reynolds number ${\textit{Re}}_{\!p}=10$ as shown in figure 2(b).

Figure 2. Effect of initial particle volume fraction on critical Rayleigh number ${\textit{Ra}}_c$ and critical wavenumber $k_c$ : (a) variation in critical Rayleigh number ${\textit{Ra}}_c$ , and (b) variation in critical wavenumber $k_c$ with initial undisturbed particle volume fraction $\varPhi _0$ for two different particle Reynolds numbers ${\textit{Re}}_{\!p}$ and particle sizes $\delta$ , and other parameters kept at $\varTheta _{\textit{pt}}=0$ , ${R}=800$ , ${E}=3385$ and ${Pr}=0.71$ for both graphs. Here, the circles represent data from Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021).

Figure 3. Effect of particle injection temperature $\varTheta _{\textit{pt}}$ on critical parameters, base-state fluid temperature and its stratification. (a) Variation in critical Rayleigh number ${\textit{Ra}}_c$ , where the dots represent the data from Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021). (b) Variation in critical wavenumber $k_c$ . (c) Variation in base-state fluid temperature. (d) Variation in unstably stratified layer thickness $\delta _{st}$ with particle injection temperature. Here, ${\textit{Re}}_{\!p}=1$ , $\delta =0.01$ , $\varPhi _0={10^{-4}}$ , ${R}=800$ , ${E}=3385$ and $ \textit{Pr}=0.71$ for all graphs.

The effect of particle injection temperature $\varTheta _{\textit{pt}}$ on the critical parameters and the base-state fluid temperature is shown in figure 3(a–d). In figure 3(a), the present data are compared with the existing work by Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021). It is observed that with an increase in $\varTheta _{\textit{pt}}$ , the critical Rayleigh number ${\textit{Ra}}_c$ increases in the initial part of figure 3(a). The explanation for this is that the particles act as the sources of heat, and by increasing their temperature, they tend to reduce the thermal stratification inside the domain, leading to an increase in the stability of the flow. However, from figure 3(c), it is clear that as the particle temperature increases beyond $\varTheta _{\textit{pt}}\gt 1$ , the strong negative base-state temperature gradients start to appear in the upper part of the domain, and favour the instability. The effect of $\varTheta _{\textit{pt}}$ on the critical wavenumber is shown in figure 3(b). As $\varTheta _{\textit{pt}}$ is increased from $-1$ to 2, $k_c$ decreases until $\varTheta _{\textit{pt}}\approx 0.5$ , then increases monotonically. This can be explained by looking at the variation in unstably stratified layer thickness $\delta _{st}$ in the base-state temperature profile shown in figure 3(c) and 3(d). For $\varTheta _{\textit{pt}}=-1$ in figure 3(c), from $z\approx -0.26$ to $-0.5$ , the fluid is unstably stratified over thickness $\delta _{st}$ , and in the remaining domain, $\varTheta _{\textit{pt}}$ has a symmetric distribution. A similar distribution exists for $\varTheta _{\textit{pt}}=2$ , but the unstably stratified layer exists near the cold top surface. The thickness of this unstably stratified layer increases for $\varTheta _{\textit{pt}}=-1$ to $\varTheta _{\textit{pt}}\approx 0$ , maintains a constant total height $\delta _{st}=1$ from $\varTheta _{\textit{pt}}\approx 0$ to $\varTheta _{\textit{pt}}\approx 1.25$ , and subsequently increases monotonically for $\varTheta _{\textit{pt}}\approx 1.25$ to $\varTheta _{\textit{pt}}=2$ . From the definition of $\delta _{st}$ , we note that the non-dimensional temperature difference across all these unstably stratified layers is 1. Hence an increase in $\delta _{st}$ value with the same temperature difference results in a lower temperature gradient, favouring stability. This explains the initial increase in ${\textit{Ra}}_c$ value in figure 3(a). From $\varTheta _{\textit{pt}}\approx 1.25$ to $\varTheta _{\textit{pt}}=2$ , the unstable stratified layer thickness reduces, maintaining the same temperature difference across its length, which increases the negative temperature gradient and favours the instability. Hence ${\textit{Ra}}_c$ reduces from $\varTheta _{\textit{pt}}\approx 1.25$ onwards. It should be noted that $k_c$ represents the length scale for the onset of convection, and it changes with particle volume fraction $\varPhi _0$ and size $\delta$ . However, in this case, we keep $\varPhi _0$ and $\delta$ constant, and vary only $\varTheta _{\textit{pt}}$ . Hence the explanation for the non-monotonic variation of $k_c$ with $\varTheta _{\textit{pt}}$ can be deduced from the variation of $\delta _{st}$ with $\varTheta _{\textit{pt}}$ shown in figure 3(d). As $\delta _{st}$ increases, the length scale at the onset of convection increases, leading to a decrease in $k_c$ , and vice versa.

4. Formulation of finite amplitude equations

In the present work, the finite amplitude expansions are given based on the analysis of Stuart (Reference Stuart1960), Yao & Rogers (Reference Yao and Rogers1992), Khandelwal & Bera (Reference Khandelwal and Bera2015), Sharma, Khandelwal & Bera (Reference Sharma, Khandelwal and Bera2018) and Aleria, Khan & Bera (Reference Aleria, Khan and Bera2024). The Fourier expansion of fluid temperature $\theta$ in separable form is

(4.1) \begin{align} \theta (x,z,t)&=\varTheta (z,\tau )\,\mathbb{E}^0+\hat {\theta }_1\left (z,\tau \right )\,\mathbb{E}^1+\hat {\theta }_2\left (z,\tau \right )\,\mathbb{E}^2+\cdots +\text{c.c.} \nonumber \\&=\mathbb{E}^0\big [\varTheta _0(z)+c_i\,|B(\tau )|^2\,\varTheta _1(z)+\mathcal{O}\big (c_i^2\big )\big ] \nonumber\\&\quad {}+\mathbb{E}^1\left [c_i^{1/2}\,B(\tau )\,\theta _{10}(z)+c_i^{3/2}\,B(\tau )\,|B(\tau )|^2\,\theta _{11}(z)+\mathcal{O}\big (c_i^{5/2}\big )\right ] \nonumber\\&\quad {}+\mathbb{E}^2\big [c_i\,B(\tau )^2\,\theta _{20}(z)+\mathcal{O}\big (c_i^2\big )\big ]+\cdots +\text{c.c.,} \end{align}

where $\mathbb{E}^j=\exp (j [\mbox{i}k_c (x-c_rt ) ])$ for $j=0,1,2$ , $k_c$ is the wavenumber corresponding to the critical Rayleigh number ${\textit{Ra}}_c$ , and $c_r$ is the real part of $c$ corresponding to the most unstable disturbance. The amplitude function $B$ will be obtained from the Landau equation, and c.c. stands for the complex conjugate of the complex-valued functions (the second and third terms).

The evolution equation for the slowly varying amplitude $B$ is derived using the method of multiple time scales. Hence two distinct time scales (a fast time scale $t$ and a slow time scale $\tau $ ) are used in the nonlinear stability analysis of the flow. As in linear stability analysis, the exponential variation of disturbance amplitude is associated with the fast time scale. As the disturbance amplitude grows exponentially and attains a finite amplitude, the nonlinear terms become important. This leads to the deviation of the temporal exponential growth of the disturbances. Hence another time scale (the slow time scale) is needed to capture the effect of nonlinearities of different orders. The slow time scale is defined as $\tau =c_it$ , with $c_i$ being the growth rate of the most unstable disturbance obtained from the linear stability analysis. Hence the derivative with respect to $t$ of a function $f(t)=F(t,\tau )$ is given by

(4.2) \begin{equation} \frac {\partial f}{\partial t}=\frac {\partial F}{\partial t}+c_i\frac {\partial F}{\partial \tau }. \end{equation}

Thus for a function of the form $F(t, \tau )$ , both $t$ and $\tau$ are treated as if they are independent variables. The justification of (4.1) is as follows. From the linear theory, it is known that the growth rate $c_i$ is linearly proportional to $ ({Ra}-{\textit{Ra}}_c )/{\textit{Ra}}_c$ , the Rayleigh number difference from the neutral curve. Using a finite amplitude stability analysis for the Rayleigh–Bénard convection (Schlüter et al. Reference Schlüter, Lortz and Busse1965) and for the non-isothermal flow in a vertical annulus, Yao & Rogers (Reference Yao and Rogers1992) showed that the square of equilibrium amplitude $A_e^2$ is proportional to $ ({Ra}-{\textit{Ra}}_c )/{\textit{Ra}}_c$ , and consequently the equilibrium amplitude is proportional to $c_i^{1/2}$ , with a proportionality constant that is a weak function of ${\textit{Ra}}$ near ${\textit{Ra}}_c$ . Hence the amplitude of the disturbance fluid temperature ${\theta }_{10}$ with wavenumber $k_c$ is of the order of $c_i^{1/2}$ . The amplitude of the disturbances of higher harmonics is chosen based on the number of interactions between lower-order harmonics. For example, the harmonic with wavenumber $2k_c$ needs two disturbances with wavenumber $k_c$ , which gives its amplitude of order $c_i$ . Similarly, the higher-order correction to the base state results from the Reynolds stress due to the interaction of $\mathbb{E}^1$ and its complex conjugate $\mathbb{E}^{-1}$ with an amplitude of order $c_i$ (Stewartson & Stuart Reference Stewartson and Stuart1971).

4.1. Derivation of the cubic Landau equation

Substitute the perturbation series similar to (4.1) for all the dependent variables in (2.2)–(2.9), and separate the harmonic components. The resulting equations are shown in § A.2. The first Landau coefficient does not depend on the higher-order harmonics ( $\mathbb{E}^3,\mathbb{E}^4$ etc.). Therefore, those terms are ignored in the series (4.1). The systems of equations at different harmonics (A9)–(A13), (A15)–(A22) and (A23)–(A30) are solved at increase orders of $c_i$ . At $\mathcal{O} (c_i^0 )$ , $\mathbb{E}^0$ harmonic equations are identical to the basic-flow equations (3.1)–(3.3). At $\mathcal{O} (c_i )^{1/2}$ , $\mathbb{E}^1$ equations are similar to that of linear stability (A1)–(A8), and the $\mathbb{E}^0$ and $\mathbb{E}^2$ do not contribute at this order. The functions $u_{x10}$ , $u_{z10}$ , $\theta _{10}$ , $v_{x10}$ , $v_{z10}$ , $\theta _{p10}$ and $\phi _{10}$ are the eigenfunctions of linear stability equations at a particular wavenumber and Rayleigh number ${\textit{Ra}}$ . At $\mathcal{O} (c_i )$ , $\mathbb{E}^0$ and $\mathbb{E}^2$ generate the system of non-homogeneous equations for basic-flow distortion functions $\varTheta _1$ , $\varTheta _{p1}$ , $V_1$ and $\varPhi _1$ , and the functions $u_{x20}$ , $u_{z20}$ , $\theta _{20}$ , $v_{x20}$ , $v_{z20}$ , $\theta _{p20}$ and $\phi _{20}$ , respectively. These equations contain the non-homogeneous part made up of known functions $u_{x10}$ , $u_{z10}$ , $\theta _{10}$ , $v_{x10}$ , $v_{z10}$ , $\theta _{p10}$ , $\phi _{10}$ and their derivatives, which are obtained at lower-order computations. At $\mathcal{O} (c_i^{3/2} )$ , the $\mathbb{E}^1$ harmonic results in a non-homogeneous system of equations with linear stability operators acting on the functions $u_{x11}$ , $u_{z11}$ , $\theta _{11}$ , $v_{x11}$ , $v_{z11}$ , $\theta _{p11}$ and $\phi _{11}$ . However, the right-hand-side non-homogeneous part contains the terms proportional to the unknown terms ${\textrm d}B/{\textrm d}\tau$ , $B$ and $B\,|B|^2$ , with coefficients depending on functions known from the lower-order analysis. The Fredholm alternative solvability condition is imposed on the non-homogeneous right-hand side of the $\mathbb{E}^1$ harmonic at $\mathcal{O} (c_i^{3/2} )$ to obtain a non-trivial solution. To impose the solvability condition, the solution to the homogeneous adjoint system (see § A.4) corresponding to the linear stability operator is required. Accordingly, the non-homogeneous right-hand side must be orthogonal to the adjoint functions $\hat {u}_x^\dagger$ , $\hat {u}_z^\dagger$ , $\hat {\theta }^\dagger$ , $\hat {v}_x^\dagger$ , $\hat {v}_z^\dagger$ , $\hat {\theta }_{p}^\dagger$ and $\hat {\phi }^\dagger$ . The orthogonality is imposed by multiplying the right-hand side of the system of linear equations of $\mathbb{E}^1$ harmonic functions at $\mathcal{O} (c_i^{3/2} )$ with $\displaystyle { (\hat {p}^\dagger , \hat {u}_{x}^\dagger ,\hat {u}_z^\dagger ,\hat {\theta }^\dagger , \,\hat {\phi }^\dagger ,\hat {v}_x^\dagger ,\hat {v}_z^\dagger ,\hat {\theta }_{\!p}^\dagger )}$ , and integrating with respect to $z$ from $-1/2$ to $1/2$ , which gives the cubic Landau equation

(4.3) \begin{equation} \frac {\textrm{d} B}{\textrm{d} \tau }=k_cB + a_1B\,|B|^2, \end{equation}

where $a_1$ is the Landau constant defined as

(4.4) \begin{equation} a_1=-\int _{-1/2}^{1/2}{\big (\mathcal{G}_x\hat {u}_{x}^\dagger +\mathcal{G}_z\hat {u}_{z}^\dagger +\mathcal{G}_\theta \hat {\theta }^\dagger +\mathcal{G}_\phi \hat {\phi }^\dagger +\mathcal{G}_{px}\hat {v}_{x}^\dagger +\mathcal{G}_{pz}\hat {v}_{z}^\dagger +\mathcal{G}_{p\theta } \hat {\theta }_{\!p}^\dagger \big )}\,{\textrm d}z, \end{equation}

which represents the first correction to the growth rate given by linear stability analysis.

By changing the variables as $A=c_i^{1/2}B$ and $\tau =c_it$ in (4.3),

(4.5) \begin{equation} \frac {\textrm{d} A}{\textrm{d} t}=k_cc_iA+a_1A^2\tilde {A}, \end{equation}

and its corresponding complex conjugate equation is

(4.6) \begin{equation} \frac {\textrm{d} \tilde {A}}{\textrm{d} t}=k_cc_i\tilde {A}+\tilde {a}_1\tilde {A}^2A, \end{equation}

where $k_cc_i$ is the growth rate from the linear stability analysis, and $A$ is the physical amplitude of the wave. Multiply (4.5) with $\tilde {A}$ and (4.6) with $A$ and add them:

(4.7) \begin{equation} \frac {\textrm{d} |A|^2}{\textrm{d} t}=2k_cc_i\,|A|^2+2\,{\textrm{Re}}\{a_1\}\,|A|^4, \end{equation}

where ${\textrm{Re}}\{a_1\}$ is the real part of $a_1$ . Equation (4.7) has an equilibrium amplitude $A_e$ as a solution if

(4.8) \begin{equation} \frac {\textrm{d} |A|}{\textrm{d} t}=0. \end{equation}

Consequently, three possible equilibrium amplitudes exist:

(4.9a,b) \begin{align} A_e=0\quad {\textrm and}\quad A_e=\pm \sqrt {-k_cc_i/{\textrm{Re}}\left \{a_1\right \}}. \end{align}

Here, $A_e=0$ represents the steady base flow, which is stable for ${\textit{Ra}}\lt {\textit{Ra}}_c$ and unstable for ${\textit{Ra}}\gt {\textit{Ra}}_c$ , and ${\textit{Ra}}_c$ is the critical Rayleigh number obtained from linear stability analysis. From (4.9), the existence of a finite amplitude solution is guaranteed if $k_cc_i$ and ${\textrm{Re}}\{a_1\}$ have opposite signs (Drazin & Reid Reference Drazin and Reid2004; Shukla & Alam Reference Shukla and Alam2011). Therefore, there are two possibilities: first, the growth rate is positive (for ${\textit{Ra}}\gt {\textit{Ra}}_c$ ), and the real part of the Landau constant is negative; second, the growth rate is negative (for ${\textit{Ra}}\lt {\textit{Ra}}_c$ ), and the real part of the Landau constant is positive. The first possibility leads to a supercritical pitchfork bifurcation, whereas the second possibility leads to a subcritical pitchfork bifurcation.

Using equilibrium amplitude definition (4.9b ), (4.7) can be rewritten as

(4.10) \begin{equation} \frac {\textrm{d} |A|^2}{\textrm{d} t}=2\,{\textrm{Re}}\left \{a_1\right \}|A|^2{\big (|A|^2-A_e^2\big )}, \end{equation}

where $A_e^2=-{ (k_cc_i )}/{\textrm{Re}}\{a_1\}$ . The solution for the above equation is

(4.11) \begin{equation} |A|^2=\frac {A_e^2}{1+\displaystyle {\left (\frac {A_e^2}{|A_0|^2}-1\right )}e^{-2k_cc_it}}, \end{equation}

where $|A_0|$ is the value of $|A|$ at $t=0$ . From (4.11), it is clear that when ${\textit{Ra}}\gt {\textit{Ra}}_c$ , as $t\to \infty$ , $|A|\to A_e$ . Hence, irrespective of the initial amplitude $|A_0|$ , the amplitude eventually tends to the equilibrium amplitude $A_e$ . Multiply (4.5) with $1/\tilde {A}$ and (4.6) with $-A/\tilde {A}^2$ , and add them:

(4.12) \begin{equation} \frac {\textrm{d} }{\textrm{d} t}{\left (\frac {A}{\tilde {A}}\right )} =2\textrm{i}A^2\, \textrm{Im} \{a_1\}, \end{equation}

where $\textrm{Im} \{a_1\}$ is the imaginary part of $a_1$ . The above equation can be rewritten as

(4.13) \begin{equation} \frac {\textrm{d}\displaystyle {\left (\frac {A}{\tilde {A}}\right )}}{\displaystyle {\left (\frac {A}{\tilde {A}}\right )}}=2{\textrm i}\,\textrm{Im} \{a_1\}\,|A|^2\,{\textrm d}t. \end{equation}

Hence the solution for $A$ is given by

(4.14) \begin{equation} \frac {A}{A_0}=\displaystyle \frac {|A|}{|A_0|}\exp \left \{\displaystyle \textrm{i}\,{\textrm{Im}} \{a_1\}\int _0^t|A|^2\,{\textrm d}t\right \}\!, \end{equation}

where $|A|{ (t )}$ is given by the square root of (4.11). The closed-form solution for amplitude function $A{ (t )}$ is obtained by integrating the above equation:

(4.15) \begin{equation} \displaystyle \frac {A(t)}{A_0}=\frac {|A|}{|A_0|}\left \{\frac {|A|}{|A_0|}\,e^{ -k_cc_it}\right \}^{{\textrm{i}\,{\textrm{Im}} \{a_1\}}/{{\textrm{Re}}\{a_1\}}}, \end{equation}

where $|A|(t)$ is obtained from (4.11).

4.2. Heat flux balance

The solution to the Landau equation given by (4.11) and (4.15) shows the existence of a steady-state solution as $t\to \infty$ . Hence in this subsection, we derive an equation for the average heat flux balance at steady state. At steady state, integrating the fluid energy (2.5) in the domain ${ (x,z )}\in \left \{{ (-0.5,0.5 )}\times { (-\pi /k_c,\pi /k_c )}\right \}$ yields

(4.16) \begin{align} &\int _{z=-0.5}^{0.5}\int _{x=-\pi /k_c}^{\pi /k_c}{\left (u_x\frac {\partial \theta }{\partial x}+u_z\frac {\partial \theta }{\partial z}\right )}\,{\textrm d}x\,{\textrm d}z \nonumber\\ &\quad =\int _{z=-0.5}^{0.5}\int _{x=-\pi /k_c}^{\pi /k_c}{\left (\frac {1}{\sqrt {{{\textit{Ra}}{\textit{Pr}}}}}{\left (\frac {\partial ^2 {\theta }}{\partial {x}^2}+\frac {\partial ^2 {\theta }}{\partial {z}^2}\right )}-\frac {{E}\phi }{{\textit{St}}_{\textit{th}}}{\left (\theta _{\!p}-\theta _{\!p}\right )}\right )}\,{\textrm d}x\,{\textrm d}z. \end{align}

Since $u_x$ , $u_z$ , $\theta$ and $\theta _{\!p}$ are periodic functions in $x$ with period $2\pi /k$ , the non-zero contribution comes only from the $\mathbb{E}^0$ harmonic. Hence (4.16) is essentially equivalent to the integration of $\mathbb{E}^0$ harmonic (A10) in $z\in { (-0.5,0.5 )}$ . Therefore, after substituting $A_e^2=c_i\,|B|^2$ , $\varTheta =\varTheta _0+A_e^2\varTheta _1$ and $\varTheta _{\!p}=\varTheta _{p0}+A_e^2\varTheta _{p1}$ in (A10), we have

(4.17) \begin{align} &\int _{z=-0.5}^{0.5}\frac {\textrm{d}^2 \varTheta }{\textrm{d}z^2}\,{\textrm d}z\\ \nonumber &={E}\sqrt {{{\textit{Ra}}{\textit{Pr}}}}\int _{z=-0.5}^{0.5}\left \{\varPhi _0\frac {\varTheta -\varTheta _{\!p}}{{\textit{St}}_{\textit{th}}} + A_e^2\varPhi _1\frac {\varTheta _0-\varTheta _{p0}}{{\textit{St}}_{\textit{th}}}+A_e^2\frac {\textrm{d} }{\textrm{d} z} {\big (\tilde {\theta }_{10}u_{z10}+\theta _{10}\tilde {u}_{z10}\big )}\right \}{\textrm d}z, \end{align}

where the integral of the last term on the right-hand side goes to zero because $u_{z10}=\theta _{10}=0$ at $z=-0.5$ and 0.5. From (A13), we have the identities

(4.18) \begin{align} \frac {\varTheta -\varTheta _{\!p}}{{\textit{St}}_{\textit{th}}}&=-v_0\frac {\textrm{d} \varTheta _{\!p}}{\textrm{d} z}+A_e^2V_1\frac {\textrm{d} \varTheta _{p0}}{\textrm{d} z}, \end{align}

(4.19) \begin{align} \frac {\varTheta _0-\varTheta _{p0}}{{\textit{St}}_{\textit{th}}}&=-v_0\frac {\textrm{d} \varTheta _{p0}}{\textrm{d} z}. \end{align}

Hence (4.17) simplifies to

(4.20) \begin{equation} [\![{\textit{Nu}}]\!]=v_0{E}\varPhi _0\sqrt {{{\textit{Ra}}{\textit{Pr}}}}\,\displaystyle [\![ \varTheta _{\!p} ]\!]+A_e^2\int _{z=-0.5}^{0.5}{\left (\varPhi _0V_1-\varPhi _1v_0\right )}\frac {\textrm{d} \varTheta _{p0}}{\textrm{d} z}\,{\textrm d}z, \end{equation}

where the bracket is defined for the function $f(z)$ as $[\![ f ]\!]=f(0.5)-f(-0.5)$ , and ${\textit{Nu}}_h={\textit{Nu}}(-0.5)=- {\textrm{d} \varTheta }/{\textrm{d} z}\big |_{z=-0.5}$ and ${\textit{Nu}}_c={\textit{Nu}}(0.5)=- {\textrm{d} \varTheta }/{\textrm{d} z}\big |_{z=0.5}$ are the Nusselt numbers (or non-dimensional heat fluxes) at the hot and cold surfaces, respectively. It can be shown that the solution to the linearised volume fraction equation is always zero (see § A.6). Hence from (A11),

(4.21) \begin{equation} \frac {\textrm{d}}{\textrm{d} z} {{\left (\varPhi _0V_1-\varPhi _1v_0\right )}}=0 \end{equation}

for $z\in { (-0.5,0.5 )}$ . Therefore, from (4.20) and (4.21), the net flux balance is given by

(4.22) \begin{equation} [\![{\textit{Nu}}]\!]=v_0{E}\varPhi _0\sqrt {{{\textit{Ra}}{\textit{Pr}}}}\, [\![ \varTheta _{\!p} ]\!], \end{equation}

where the proportionality constant $v_0{E}\varPhi _0\sqrt {{{\textit{Ra}}{\textit{Pr}}}}$ is equal to ${L}^2/l$ . Here, $L$ and $l$ are the non-dimensional length scales given by (3.11) and (3.12), respectively. Thus the alternate form of the above equation is given by

(4.23) \begin{equation} [\![{\textit{Nu}}]\!]=\displaystyle [\![ Q_{\!p}^{\prime \prime } ]\!], \end{equation}

where $Q_{ph}^{\prime \prime }=Q_{\!p}^{\prime \prime }(-0.5)=(l/{L}^2)\,\varTheta _{\!p}(-0.5)$ and $Q_{pc}^{\prime \prime }=Q_{\!p}^{\prime \prime }(0.5)=(l/{L}^2)\,\varTheta _{\!p}(0.5)$ are the non-dimensional particle sensible heat fluxes at hot and cold surfaces, respectively. The physical significance of (4.23) is that at steady state, the net heat flux released from the hot and cold surfaces is equal to the net sensible heat exchanged by the particles. Therefore, for the single-phase and particle-laden flows without thermal energy two-way coupling, the heat flux at the hot surface is exactly the heat flux at the cold surface. However, for the problems with thermal energy two-way coupling, the two phases exchange the heat by following (4.23), and heat fluxes $Nu$ at hot and cold need not be equal.