1. Introduction
Thermally driven convection is of critical importance to a wide variety of natural phenomena and engineering applications (Busse Reference Busse1978; Cardin & Olson Reference Cardin and Olson1994; Rahmstorf Reference Rahmstorf2000; Hartmann, Moy & Fu Reference Hartmann, Moy and Fu2001). In many practical scenarios, the emergence of thermal convection is usually accompanied by dispersed particulate phases, such as volcanic ash clouds, dusty pollutants spreading into the atmosphere (Shaw Reference Shaw2003; Lube et al. Reference Lube, Breard, Esposti-Ongaro, Dufek and Brand2020) and the solidified crystals nucleated in magma oceans (Patočka et al. Reference Patočka, Tosi and Calzavarini2022). In particle-laden thermal convection, there exists a special scenario where particles are not only passively advected by the flow, but they also actively drive the convection through heat release. A canonical example is the flow in the particle solar receiver, where the working fluid is laden with small dispersed particles to absorb the concentrated solar radiation (Pouransari & Mani Reference Pouransari and Mani2018; Banko et al. Reference Banko, Villafañe, Kim and Eaton2020). In such flows, dispersed inertial particles harvest energy from external radiation, acting as moving heat sources and continuously heating the surrounding fluid, thus they can affect the flow dynamics through both thermal and momentum interactions between particles and turbulence. These interphase interactions make the thermal convection driven by heat-releasing particles (particle-heating convection (PHC)) much more sophisticated and therefore poorly understood hitherto, compared with its single-phase counterparts driven purely by internal heat sources or wall temperature differences (Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000; Lohse & Xia Reference Lohse and Xia2010).
In contrast to the single-phase situations, the particle-laden thermal convection is known to be signatured by intricate particle–fluid interactions. Extensive studies have investigated the convection dynamics modified by inertial particles that involve only the momentum coupling, mainly in the classical Rayleigh–Bénard convection (RBC) driven by the temperature difference between bottom and top walls (Zonta, Marchioli & Soldati Reference Zonta, Marchioli and Soldati2008; Ling, Balachandar & Parmar Reference Ling, Balachandar and Parmar2016; Pan et al. Reference Pan, Dong, Zhou and Shen2022; Xu, Xu & Xi Reference Xu, Xu and Xi2024). Nevertheless, striking modifications have been observed for the RBC laden with isothermal particles. Specifically, Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021) found that particles introduced at the top plate have a stabilizing effect on RBC, with the stability threshold magnitude increasing up to several orders. In the RBC laden with isothermal inertial particles in the range of
$10^6 \leqslant Ra \leqslant 10^8$
(Yang et al. Reference Yang, Zhang, Wang, Dong and Zhou2022b
), where
${\textit{Ra}}$
is the Rayleigh number measuring the strength of buoyancy, the global heat transfer and turbulent momentum transfer are remarkably enhanced for a Stokes number
${\textit{St}}_p=8\times 10^{-3}$
. Through interphase coupling, the flow carries particles into the near-wall region and consequently, particles preferentially concentrate in the streaming zones with downward fluid. For a weakly turbulent RBC of
${\textit{Ra}}=2\times 10^6$
, Oresta & Prosperetti (Reference Oresta and Prosperetti2013) reported strong effects of the two-way momentum coupling on heat transfer, characterized by a notable increase of Nusselt number
$({\textit{Nu}})$
versus the increasing size of point-like settling particles. The mechanism is attributed to the preferential concentration of falling particles in the sinking stream, which drags and accelerates the local fluid downward, thereby promoting heat transfer. Via direct numerical simulations (DNS), Alards et al. (Reference Alards, Kunnen, Clercx and Toschi2019) have studied the dynamics of thermally expandable point particles in RBC for
${\textit{Ra}}=2\times 10^7$
, where the particles are considered to have a thermal expansion coefficient larger than that of the background fluid. It is found that the particles with a larger thermal expansion coefficient and inertia would accumulate in the centres of Rayleigh–Bénard cells. A similar problem was experimentally investigated by Hu et al. (Reference Hu, Wang, Jia, Zhong and Zhang2021) for rod-like particles with
${\textit{Ra}}$
ranging from
$3\times 10^8 - 4\times 10^9$
. These particles act as ‘active’ mixers of the velocity and temperature fields across the boundary layers, enhancing heat transport at relatively higher
${\textit{Ra}}$
. For a conceptualized ‘active particle’ turbulence, i.e. a water-based turbulent convection system with a minute concentration of hydrofluoroether introduced, the biphasic dynamics is found to enhance the global heat transfer efficiency up to
$500 \, \%$
(Wang et al. Reference Wang, Mathai and Sun2019, Reference Wang, Mathai and Sun2020b
).
In contrast to isothermal particles, heat-releasing particles play a distinctly important role in thermal convection, as thermal coupling is introduced and comes into play in the particle–fluid interactions. That is, they can serve not only as moving heat sources to the surrounding fluid but also as inducers of intense local driving forces (buoyancy). These effects become much more profound when preferential concentration occurs in various flow systems. Specifically, the preferential concentration of heat sources significantly alters turbulent flow dynamics and interphase heat transfer efficiency (Kuerten, Van der Geld & Geurts Reference Kuerten, Van der Geld and Geurts2011; Pouransari & Mani Reference Pouransari and Mani2018; Rahmani et al. Reference Rahmani, Geraci, Iaccarino and Mani2018). For the PHC in a cubic domain imposed with periodic boundary conditions, the pioneering work of Zamansky et al. (Reference Zamansky, Coletti, Massot and Mani2014) uncovered the feedback loop mechanism resulting from the multiway interphase interaction, showing that the fluid temperature fluctuations are intensified due to the preferential concentration. In their following work (Zamansky et al. Reference Zamansky, Coletti, Massot and Mani2016), this feedback loop is classified into two different regimes depending on the particle inertia. For very small particle inertia, many microscale thermal plumes are generated independently around particles and give rise to the macroscopic flow behaviour. For large particle inertia, particle clusters take form and their dynamics are influenced by both larger structures and small-scale dynamics. Importantly, these clusters provide a feedback modification on the large-scale flow dynamics, via the rearrangement of the multiscale distribution of heat sources. Similar modification of large turbulence structures is also observed in forced homogeneous turbulence (Frankel et al. Reference Frankel, Pouransari, Coletti and Mani2016), where the particle clusters modify the flow dynamics by shedding hot buoyant plumes. Meanwhile, the mean settling velocity is significantly reduced. In summary, preferential concentration of heat-releasing particles has profound effects on large-scale flow structures.
As is well known, the large-scale circulation (LSC) is regarded as the most important large-scale flow structure in thermal convection systems, identified as well-shaped flow patterns in the background of turbulent fluctuations (Brown & Ahlers Reference Brown and Ahlers2007; Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Xia Reference Xia2011; Chong et al. Reference Chong, Wagner, Kaczorowski, Shishkina and Xia2018; Xu et al. Reference Xu, Chen, Wang and Xi2020). Most existing research on the LSC dynamics is mainly based on the classical single-phase RBC systems (Bodenschatz et al. Reference Bodenschatz, Pesch and Ahlers2000; Lohse & Xia Reference Lohse and Xia2010; Ecke & Shishkina Reference Ecke and Shishkina2023), and numerous studies have highlighted the pivotal role of the LSC in dominating the heat and mass transfer. For example, introducing weak perturbations to the LSC can significantly modify the global transfer efficiency and turbulence intensity (Wang et al. Reference Wang, Xu, Xia, Wan and Sun2017, Reference Wang, Wan, Yan and Sun2018a
, Reference Wang, Chen, Xu and Xi2024b
). On the other hand, the topology and dynamics of the LSC are highly sensitive to key controlling parameters (e.g. Rayleigh number
${\textit{Ra}}$
and Prandtl number
${\textit{Pr}}$
) and flow configurations (e.g. the geometrical set-up and boundary conditions) (Wang, Mathai & Sun Reference Wang, Mathai and Sun2020b
). This has been demonstrated in a canonical cylindrical RBC cell of aspect ratio 1, by which rich dynamic features are observed including azimuthal rotations, occasional cessations, torsional oscillations and reversals (Sreenivasan, Bershadskii & Niemela Reference Sreenivasan, Bershadskii and Niemela2002; Xi, Zhou & Xia Reference Xi, Zhou and Xia2006; Funfschilling, Brown & Ahlers Reference Funfschilling, Brown and Ahlers2008; Mamykin et al. Reference Mamykin, Kolesnichenko, Pavlinov and Khalilov2018; Wei Reference Wei2021; Chen et al. Reference Chen, Xu, Xia and Xi2023). These dynamical features, especially the LSC reversal, provide insightful perspectives into the complex convection system.
To date, limited research attempts have been directed to examining the influence of heat-releasing particles on the LSC dynamics, the commensurate heat transfer and the feedback particle behaviours, although changes in heat source distribution are known to cause intriguingly different flow physics in thermal convection (Tasaka & Takeda Reference Tasaka and Takeda2005; Kuznetsov & Nield Reference Kuznetsov and Nield2016; Bouillaut et al. Reference Bouillaut, Lepot, Aumaître and Gallet2019). This may be attributed to the fact that the LSC appearing in the classic RBC of the enclosed domain is absent when thermal convection driven by heated inertial particles is considered in a periodic cubic domain (Zamansky et al. Reference Zamansky, Coletti, Massot and Mani2014, Reference Zamansky, Coletti, Massot and Mani2016). To this point, by applying non-slip boundary conditions for the top and bottom walls and considering the horizontal directions to be periodic, Du & Yang (Reference Du and Yang2022) have observed the LSC manifested as convection rolls in the PHC. However, they did not conduct a comprehensive study regarding the underlying mechanisms of the LSC dynamics, while focused mainly on the dependence of response parameters on
${\textit{Ra}}$
and
${\textit{St}}_p$
. Of interest, Du & Yang (Reference Du and Yang2022) also identified a regime with most of the heat-releasing particles accumulating into the top boundary layer at intermediate Stokes number
${\textit{St}}_p \! \sim \! O(1)$
.
Considering the above facts, the present study is dedicated to a careful exploration of the novel LSC dynamics, as well as the underlying mechanism of the interphase momentum and thermal coupling in the PHC. To this end, extensive DNS are designed and performed in a cubic domain confined by solid walls. The simulations consider wide ranges of particle inertia (manifested here by the density ratio between particle and fluid) and thermal buoyancy strength (related to the heating power exerted on particles). To reliably resolve the effect of preferential concentration on the LSC dynamics, the well-established four-way coupling method is implemented (Squires & Eaton Reference Squires and Eaton1991; Eaton & Fessler Reference Eaton and Fessler1994; Tom & Bragg Reference Tom and Bragg2019; Liao et al. Reference Liao, Chen, Wan, Liu and Lu2024), taking into account both momentum and heat transfer between particles during the collision process. The intriguing correlations underlying the LSC regime transition, the commensurate flow statistics and particle distribution features are discussed. Furthermore, a mechanistical interpretation is proposed for the onset of the LSC regime transition, based on a quantitative description of particle centrifugal motion.
This paper is organized as follows. In § 2, the problem set-up, governing equations, dimensionless parameters and numerical methods are described. Results and discussion are detailed in § 3. Finally, concluding remarks are summarized in § 4.
Schematic of thermal convection driven by heat-releasing point particles in a cube.
2. Mathematical models and numerical procedure
2.1. Governing equations for fluid phase
As sketched in figure 1, thermal convection driven by heat-releasing particles (coloured in red) is confined in a cube of edge length
$L$
. The heat-releasing particles of density
$\rho _p$
have a total population
$N_p$
and each particle is subjected to a given external heating power
$Q_p^*$
. In other words, a total external heating power
$N_p Q_p^*$
is introduced into this two-phase flow system via heating individual particles radiatively, which in the meantime conduct heat transfer with the surrounding fluid (Yang et al. Reference Yang, Wan, Zhou and Dong2022a
). The top and bottom walls (coloured in azure) are imposed with no-slip and isothermal boundary conditions for the fluid phase, and the four sidewalls (coloured in grey) with no-slip and adiabatic boundary conditions. The fluid carrying the dispersed particles takes the Oberbeck–Boussinesq approximation, with density
$\rho$
holding a linear dependence on temperature
$T$
, namely
$\rho =\rho _f[1-\beta (T-T_0)]$
, where
$\beta$
and
$\rho _f$
are the thermal expansion coefficient and density of fluid at a reference temperature
$T_0$
, respectively. The non-dimensional governing equations for the fluid phase read as
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}& = 0, \end{align}
\begin{align} \partial _t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} + \boldsymbol{\nabla }p &= \sqrt {\frac {\textit{Pr}}{\textit{Rr}}} {\nabla} ^2 \boldsymbol{u} + T \boldsymbol{e}_z + \boldsymbol{F}_{\kern-1.5pt f} , \end{align}
\begin{align} \partial _t T+ \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }T &= \sqrt {\frac {1}{Pr Rr}} {\nabla} ^2 T + Q_{\kern-1.5pt f} , \end{align}
where
$\boldsymbol{u}=(u_x, u_y, u_z)$
and
$p$
denote the velocity and pressure of the carrier phase, respectively;
$\boldsymbol{F}_{\kern-1.5pt f}$
and
$Q_{\kern-1.5pt f}$
represent the interphase momentum and heat exchange rate, respectively. More details of
$\boldsymbol{F}_{\kern-1.5pt f}$
and
$Q_{\kern-1.5pt f}$
can be found in Appendix A;
$\boldsymbol{e}_z$
is the unit vector of the vertical (
$z$
) direction, opposite to the gravitational acceleration of magnitude
$g$
in the Eulerian coordinates
$\boldsymbol{x}=(x,y,z)$
. The key controlling parameters include the Prandtl number (
$Pr$
) and the Rayleigh–Robert number (
${\textit{Rr}}$
), defined, respectively, as
\begin{equation} {\textit{Pr}}=\frac {\nu }{\kappa } \quad \textrm {and} \quad {\textit{Rr}}=\frac {g\beta N_p Q_p^* L^2}{\nu \kappa ^2 \rho _f C_{\textit{pf}}}, \end{equation}
where
$\nu$
,
$\kappa$
and
$C_{\textit{pf}}$
are the fluid kinematic viscosity, thermal diffusivity and specific heat capacity, respectively. Different from the
${\textit{Ra}}$
used in the classical RBC, here the driving force is measured by
${\textit{Rr}}$
, which quantifies the total external heating power (
$N_p Q_p^*$
) input via the particle phase. A similar definition is also used in the single-phase internally heated convection (IHC) (Goluskin Reference Goluskin2016; Kuznetsov & Nield Reference Kuznetsov and Nield2016; Wang, Lohse & Shishkina Reference Wang, Lohse and Shishkina2021). It should be noted that the IHC is characterized by fluid temperature that is always higher in the bulk than the wall temperature, resulting in a stable and an unstable temperature stratification at the bottom and top walls, respectively. As a consequence, the buoyancy force acts to break the up–down symmetry by stabilizing the flow near the bottom wall while driving the flow near the top one. Hereafter, the fluid density
$\rho _f$
, the edge length
$L$
and the free-fall velocity
$U =\sqrt {g\beta T_c L}$
are used as the characteristic scales (Xu et al. Reference Xu, Xu, Jiang and Xi2022; Wondrak et al. Reference Wondrak, Sieger, Mitra, Schindler, Stefani, Vogt and Eckert2023). The characteristic temperature (
$T_c$
) is obtained as (Goluskin Reference Goluskin2016; Wang et al. Reference Wang, Lohse and Shishkina2021)
\begin{equation} T_c = \dfrac {N_p Q^*_p}{\kappa \rho _f C_{\textit{pf}}L}. \end{equation}
2.2. Governing equations for particle phase
The dispersed particles are considered to have identical physical properties and a sufficiently small diameter
$d_p$
, and thus can be regarded as point heat sources. They are tracked from the Lagrangian perspective with four-way coupling realized in the context of the discrete element method (DEM) (Schmeeckle Reference Schmeeckle2014). The particle motion and thermal dynamics are governed by the following non-dimensional equations:
\begin{align} \frac {\textrm {d} \boldsymbol{v}_i}{\textrm {d}t} &= \sum _j \boldsymbol{F}_{\textit{ij}}^c + \frac {C_d {\textit{Re}}_p}{St_p}(\boldsymbol{u}_i -\boldsymbol{v}_i) + \frac {1}{Fr^2} \boldsymbol{e}_z, \end{align}
\begin{align} \frac {\textrm {d} T_i}{\textrm {d}t} &= \sum _j Q_{\textit{ij}}^c + \frac {{\textit{Nu}}_p}{St_{\theta }}(T_{\textit{fi}} -T_i) +Q_p, \end{align}
where
$\boldsymbol{v}_i = (v_x, v_y, v_z)_i$
is the velocity and
$T_i$
the temperature of the
$i$
th particle;
$\boldsymbol{F}^{c}_{\textit{ij}}$
and
$Q^{c}_{\textit{ij}}$
are, respectively, the collision force and heat exchange rate between the
$i$
th and
$j$
th particles, which are modelled by the soft-sphere approach (Brilliantov et al. Reference Brilliantov, Spahn, Hertzsch and Pöschel1996; Li & Mason Reference Li and Mason2000; Bodhisattwa, Fernando & Silvina Reference Bodhisattwa, Fernando and Silvina2006; Ai et al. Reference Ai, Chen, Rotter and Ooi2011), see more details in Appendix A;
$T_{\textit{fi}}$
is the fluid temperature at the position of particle
$i$
. Note here,
$Q_p$
is the dimensionless heating power prescribed on each particle and thus the real control parameter varying with
${\textit{Rr}}$
. The particle–wall interactions employ the same model as that for particle–particle collision, so that all particles are confined within the cube. The drag coefficient
$C_d$
is calculated by (Zhu et al. Reference Zhu, Zhou, Yang and Yu2007a
)
\begin{equation} \begin{aligned} C_d= \begin{cases} 24\big(1+0.15{\textit{Re}}_p^{0.687}\big)/{\textit{Re}}_p,& \,{\textit{Re}}_p\lt 1000, \\ 0.44, & \,{\textit{Re}}_p \geqslant 1000. \end{cases} \end{aligned} \end{equation}
The particle Reynolds and Nusselt numbers, i.e.
${\textit{Re}}_p$
and
${\textit{Nu}}_p$
, are defined as
\begin{equation} {\textit{Re}}_p = \frac {d_p \left | \boldsymbol{u}_i -\boldsymbol{v}_i \right |}{\nu } \quad \textrm {and} \quad {\textit{Nu}}_p=\frac {h_{\textit{fi}} d_p}{\kappa \rho _f C_{\textit{pf}}}, \end{equation}
respectively, where
$h_{\textit{fi}}$
denotes the heat transfer coefficient between particle
$i$
and the surrounding fluid. Here
${\textit{Nu}}_p$
is evaluated by the widely used model developed by Ranz & Marshall (Reference Ranz and Marshall1952a
,
Reference Ranz and Marshallb
),
${\textit{Nu}}_p= 2+0.6{\textit{Re}}_p^{0.5}{\textit{Pr}}^{0.33}$
. The key parameters that commonly quantify the momentum and thermal inertia effects of particles are the momentum (
${\textit{St}}_p$
) and thermal (
$St_{\theta }$
) Stokes numbers, respectively. They are defined as
\begin{equation} St_p = \frac {\tau _p}{L/U} \quad \textrm {and} \quad St_{\theta } = \frac {\tau _{\theta }}{L/U}, \end{equation}
based on the momentum (
$\tau _p$
) and thermal (
$\tau _\theta$
) relaxation times
\begin{equation} \tau _p=\frac {\hat {\rho }_r d_p^2}{18 \nu } \quad \textrm {and} \quad \tau _{\theta } = \frac {d_p^2 \hat {\rho }_r C_{\textit{pp}}}{\kappa C_{\textit{pf}}}, \end{equation}
with
$C_{\textit{pp}}$
being the specific heat of particles. Inherently, both
${\textit{St}}_p$
and
$St_{\theta }$
linearly depend on the density ratio
$\hat {\rho }_r$
, which is defined as
\begin{equation} \hat {\rho }_r = \rho _p/\rho _f. \end{equation}
The strength of the gravitational settling effect is characterized by the Froude number
$Fr$
, which is an independent parameter defined as
\begin{equation} Fr = \sqrt {\dfrac {\beta \hat {\rho }_r N_p Q_p^*}{(\hat {\rho }_r-1)\kappa \rho _f C_{\textit{pf}}L}}. \end{equation}
According to Zamansky et al. (Reference Zamansky, Coletti, Massot and Mani2016) and Du & Yang (Reference Du and Yang2023), when
$Fr^2 \gg St_p$
, heat-releasing particles can be advected in the entire domain without notable settling. It is therefore common to neglect the gravitational settling effect of point particles at high
$Fr$
in prior studies (Zamansky et al. Reference Zamansky, Coletti, Massot and Mani2014; Lee & Lee Reference Lee and Lee2015; Du & Yang Reference Du and Yang2022; Liao et al. Reference Liao, Chen, Wan, Liu and Lu2024). Therefore, we consider only the simulations in which
$Fr$
is assumed sufficiently large so that the gravitational settling effect is ignored.
2.3. Parameter settings and numerical details
Special interest is directed to explore the intriguing influence of thermal buoyancy and particle inertia on the LSC pattern and the related flow responses, in the PHC within a cube containing
$N_p=10^5$
particles. To this end, the DEM-based DNS is performed for
${\textit{Rr}}$
and
$\hat {\rho }_r$
over three orders, specifically
${\textit{Rr}}=4.97\times 10^{5}= 4.97 \times 10^{8}$
and
$\hat {\rho }_r=1- 1000$
, and
${\textit{Pr}}=1$
. Other parameter settings, including particle Stokes numbers
${\textit{St}}_p$
and
$St_\theta$
are listed in table 1 of Appendix A and table 2 of Appendix C, which also summarize the dimensional and dimensionless parameters involved.
Following Metzger et al. (Reference Metzger, Strässle, Girardin, Conzelmann and Müller2022) and Xie et al. (Reference Xie, Hu, Zhu, Yu and Pähtz2023a
,
Reference Xie, Zhu, Hu, Yu and Panb
), the DNS-DEM approach is realized by use of CFDEMcoupling, an open source code developed and well validated by Kloss et al. (Reference Kloss, Goniva, Hager, Amberger and Pirker2012) and Hager et al. (Reference Hager, Kloss, Pirker and Goniva2014). In the CFDEMcoupling framework, the fluid phase is solved using OpenFOAM (Jasak, Jemcov & Tukovic Reference Jasak, Jemcov and Tukovic2007; Jasak Reference Jasak2009), and the particle phase using LIGGGHTS (Mubin, Li & Plimpton Reference Mubin, Li and Plimpton2021). Specifically, for the fluid phase, spatial discretization is conducted by a second-order central difference scheme, and time integration by the second-order implicit Euler scheme. The fluid velocity and pressure are solved using the PIMPLE algorithm (Jasak Reference Jasak2009). For the particle phase, time integration is accomplished by a well-established second-order scheme, the velocity Verlet algorithm (Martys & Mountain Reference Martys and Mountain1999). According to Liao et al. (Reference Liao, Chen, Wan, Liu and Lu2024), the second-order Lagrangian interpolation is employed to obtain the variables at the particle position using surrounding Eulerian grids, and the source terms
$\boldsymbol{F}_{\kern-1.5pt f}$
and
$Q_{\kern-1.5pt f}$
are smoothed and deposited towards the eight neighbouring Eulerian grids around each particle.
All the present calculations begin with the corresponding fully developed states obtained for the IHC at the same
${\textit{Rr}}$
and
${\textit{Pr}}$
values. Heat-releasing particles are introduced randomly into the flow fields, with velocity and temperature initialized using the surrounding fluid variables. A stretched mesh is employed near the walls to fully resolve the essential turbulence dynamics. For
${\textit{Rr}}\leqslant 1.49\times 10^{8}$
, the mesh size used is
$N_x\times N_y \times N_z=120\times 120\times 160$
; while for
${\textit{Rr}}=4.97\times 10^{8}$
, a finer one of
$160\times 160\times 384$
is adopted to ensure grid resolution. The adequacy of grid resolution is confirmed by evaluating the Kolmogorov (
$\eta _u$
) and Batchelor (
$\eta _T$
) scales, i.e.
$\eta _u ={\textit{Pr}} ^{1/2} (Rr\langle u_z T \rangle )^{-1/4}$
and
$\eta _T = \eta _u {\textit{Pr}} ^{-1/2}$
(Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010; Goluskin & Van der Poel Reference Goluskin and Van der Poel2016; Du & Yang Reference Du and Yang2022). We obtain
$\eta _u, \eta _T\gt 0.01$
for all the cases, and this assures that grid resolution is sufficiently fine, i.e.
$\text{max} \{\unicode{x1D6E5} x, \unicode{x1D6E5} y, \unicode{x1D6E5} z \}\leqslant 0.008$
(Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010). More numerical details are listed in table 2 of Appendix C, and a grid convergence check is provided in Appendix B. Time integration is performed with a small time step of
$\unicode{x1D6E5} t_f \sim 10^{-2}$
for the fluid phase, and a much smaller one of
$\unicode{x1D6E5} t_p = \unicode{x1D6E5} t_f/2000 - \unicode{x1D6E5} t_f/50$
for the particle phase. In that, we have
$\unicode{x1D6E5} t_p \lt \tau _{p}, \tau _{\theta }$
to reliably resolving the fluid–particle interaction (Li, Xu & Thornton Reference Li, Xu and Thornton2005; Schmeeckle Reference Schmeeckle2014; Xie et al. Reference Xie, Zhu, Hu, Yu and Pan2023b
). Long-time calculations of typically
$t \geqslant 2200$
are conducted to achieve the stationarily steady states, and after that, turbulence statistics are calculated by ensemble averaging over 2000 snapshots of flow fields.
Numerical analyses are performed based on key statistics averaged in various contexts. Based on the time averaging (
$\langle \, \boldsymbol{\cdot }\,\rangle _{t}$
) that is conducted over a time interval
$\unicode{x1D6E5} t$
, the horizontal-time (HT) averaging (
$\,\overline {\,\boldsymbol{\cdot }\,}\,$
) is conducted further over a given
$(x,y)$
plane at position
$z$
, and the volume-time (VT) averaging (
$\langle \boldsymbol{\cdot }\rangle _\varPsi$
) is conducted further over a volume
$V_\varPsi$
that meets a given condition
$\varPsi$
. For example,
$\langle \boldsymbol{\cdot }\rangle _{z\gt 0.9}$
represents that the VT averaging is performed in the domain of
$z\gt 0.9$
. As two important examples, the VT averaging over the entire cube is denoted as
$\langle \, \boldsymbol{\cdot }\,\rangle$
, and the averaging over all the particle cells (defined below in § 3.3) that have the same minimum distance
$d(\equiv \text{min} \{ x, \, y, \, 1-y, \, 1-y \})$
to any one of the four sidewalls is denoted by
$\widetilde {\boldsymbol{\cdot }}$
. For a given variable
$f(\boldsymbol{x},\tau )$
, the above averaging operations are defined as
\begin{align} \langle f \rangle _{t} &= \frac {1}{\unicode{x1D6E5} t}\int _{\unicode{x1D6E5} t} f(\boldsymbol{x},\tau ) \textrm {d} \tau , \end{align}
\begin{align} \overline {f}(z)&=\frac {1}{L^2} \iint \langle f \rangle _{t} \, \textrm {d} x \textrm {d} y, \end{align}
\begin{align} \langle f \rangle _{\varPsi } &= \frac {1}{V_{\varPsi }} \iiint \langle f \rangle _{t} \, \textrm {d} x \textrm {d} y\textrm {d} z, \end{align}
\begin{align} \widetilde {f}(d) &= \frac {1}{V_{\varPi }} \int _{\varPi } \langle f \rangle _{t} \, \textrm {d} \varPi , \end{align}
where
$V_{\varPi }$
denotes the total volume of all the particle cells that have the same
$d$
value.
3. Results and discussion
3.1. The LSC regime transition
Three typical LSC regimes have been identified based on the time-averaged velocity fields, namely, the three-dimensional streamlines and the contour plots of vertical velocity
$\langle u_z \rangle _t$
, see figure 2 for representative visualization at
${\textit{Rr}}=4.97\times 10^7$
. Specifically, as viewed in figure 2(
$a$
,
$e$
), a bulk-flow-up (BFU) regime is observed at
$\hat {\rho }_r=1$
. This flow regime is characterized by the well-structured LSC that dominates the entire domain, with fluid in the bulk rising up and in the four corners sinking down. Therefore, it is manifested as two dominant, regular and nearly symmetric counter-rotating roll cells in the vertical diagonal plane.
Continuous increase of
$\hat {\rho }_r$
leads to an intriguing flow regime transition, demonstrating the modulating effects of particle inertia on flow topology selection. As shown in figure 2(
$b$
), the LSC at
$\hat {\rho }_r=10$
consists of basically three distorted roll cells in the vertical diagonal plane: a very large one occupying almost the entire domain and two much smaller ones residing in the corners. The dominant large roll cell has upflow and downflow in the opposite diagonal corners, similar to the observation of Qiu & Tong (Reference Qiu and Tong2001). Consequently, when viewed in the midheight plane, the flow is patterned as a bulk composed of both rising-up and sinking-down fluids, with upflow concentrated in one corner and downflow in the other three corners (see figure 2
$f$
). Therefore, it is classified as the ‘marginal regime’ hereafter to summarize the above essential flow features.
Of interest, the LSC transitions to the bulk-flow-down (BFD) regime as
$\hat {\rho }_r \geqslant 50$
. As illustrated in figure 2(
$c,d$
), the LSC in the BFD regime is quite similar to that in the BFU regime, but with sinking-down fluid occupying the bulk region. In other words, the flow regime transition demonstrates a particle-inertia-induced LSC reversal. The flow topology is reorganized as only two regular counter-rotating roll cells in the vertical diagonal plane. Different from the BFU regime, the BFD regime has rising-up flow not only in the four corners but also near the sidewalls (see figure 2
$g,h$
). One should keep in mind that this LSC reversal arises due to the increasing inertia (
$\hat {\rho }_r$
here) of heat-releasing particles, namely, a flow regime transition that happens as a statistical (time-averaged) behaviour in the parameter plane of (
${\textit{Rr}}, \hat {\rho }_r$
). To this end, it is essentially different from the LSC reversal one might recall for the classic RBC problems. Specifically, the LSC reversal extensively studied in the RBC is inherent to the nonlinear flow dynamics and manifests as a transient switch of instantaneous flow patterns commonly at fixed parameter settings (Sugiyama et al. Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010; Ni, Huang & Xia Reference Ni, Huang and Xia2015; Wang et al. Reference Wang, Lai, Song and Tong2018b
).
Typical flow regimes for
${\textit{Rr}}=4.97 \times 10^7$
:
$(a,e)$
the BFU regime,
$(b,f)$
the marginal regime and
$(c,d,g,h)$
the BFD regime, which are identified based on the time-averaged velocity fields illustrated by
$(a,b,c,d)$
the streamlines and
$(e,f,g,h)$
the contours of the vertical velocity
$\langle u_z\rangle _t$
at
$y=0.5$
. For a clear illustration, the streamlines are shown in a vertical diagonal section and coloured by
$\langle u_z\rangle _t$
.
To further characterize the flow regime transition, two first-order central moments of time-averaged vertical velocity, i.e.
$M_x$
and
$M_y$
are introduced as quantitative indicators. This is inspired by the idea well-established to effectively diagnose the LSC reversal by the sign change of global angular momentum (Sugiyama et al. Reference Sugiyama, Calzavarini, Grossmann and Lohse2009, Reference Sugiyama, Ni, Stevens, Chan, Zhou, Xi, Sun, Grossmann, Xia and Lohse2010). Specifically,
$M_x$
and
$M_y$
are defined in the midheight plane (
$y=0.5$
) as
\begin{align} M_x &= \iint \frac {\langle u_z \rangle _t}{|\langle u_z \rangle _t|_{max}} |x-0.5| \, \textrm {d} x \textrm {d} y, \\[-12pt] \nonumber \end{align}
\begin{align} M_y &= \iint \frac {\langle u_z \rangle _t}{|{\langle u_z \rangle _t}|_{max}} |y-0.5| \, \textrm {d} x \textrm {d} y. \end{align}
Clearly,
$M_x$
and
$M_y$
obtain the limiting values of
$\unicode{x1D6E5} = 0.25$
and
$-0.25$
when
$\langle u_z \rangle _t$
is positively and negatively uniform in the midheight plane, respectively. For the present study, a threshold value of
$\epsilon =0.08|\unicode{x1D6E5} |=0.02$
is chosen to identify the flow regimes. As depicted in figure 3, the BFU regime is identified as
$M_x, {M_y}\lt -\epsilon$
while the BFD regime satisfies
$M_x, {M_y}\gt \epsilon$
, and all other situations fall into the marginal regime. Based on this, the LSC regime transition observed above is easily diagnosed in the (
${\textit{Rr}},\, \hat {\rho }_r$
) plane by the variations of (
$M_x,\, {M_y}$
). Interestingly, it is indicated in figure 3
$(a,\, b)$
that the values of (
$M_x,\, {M_y}$
) in the marginal regime are sensitive to the change of
$\hat {\rho }_r$
. This confirms that the flow regime transition is driven basically by increasing particle inertia.
Variations of
$(a)$
$M_x$
and
$(b)$
$M_y$
versus density ratio
$\hat {\rho }_r$
at different
${\textit{Rr}}$
. Accordingly, the parameter settings explored are classified by
$M_x, {M_y} \lt -\epsilon$
for the BFU regime,
$M_x, {M_y} \gt \epsilon$
for the BFD regime and other situations for the marginal regime.
Phase diagram in the parameter space of (
${\textit{Rr}},\hat {\rho }_r$
) with flow regimes classified by the (
$M_x$
,
$M_y$
) values. The BFU, BFD and marginal regimes are denoted in sequence by the red squares, the blue circles and the green triangles. The black dashed line indicates a power fitting scaled as
$-2/5$
.
Based on the variations of (
$M_x, \, {M_y}$
), a phase diagram of typical flow regimes is obtained in the parameter space of (
${\textit{Rr}},\, \hat {\rho }_r$
), illustrating clearly the flow regime transition (see figure 4). It is indicated that the particle inertia (
$\hat {\rho }_r$
), rather than the thermal buoyancy (
${\textit{Rr}}$
), is more critical in triggering the LSC regime transition. In particular, for all
${\textit{Rr}}$
, gradually increasing
$\hat {\rho }_r$
leads to the transition from the BFU regime to the marginal regime, and eventually to the BFD regime. However, by increasing
${\textit{Rr}}$
at a given
$\hat {\rho }_r$
, such a transition may not always happen. For example, the flow maintains the BFU and BFD regimes at
$\hat {\rho }_r=1$
and 1000, respectively, despite the increase of
${\textit{Rr}}$
up to three orders of magnitude. Therefore, the LSC regime transition is absent for heat-releasing particles of the smallest or largest inertia considered. Intriguingly, the marginal regimes indicate a power fitting scaled as
$-2/5$
in the (
${\textit{Rr}},\, \hat {\rho }_r$
) plane, for which the underlying mechanism will be discussed in § 3.4.
3.2. Commensurate flow statistics
It is well known that the LSC is dominantly important to the efficient heat and momentum transfer in thermal convection systems (Stevens et al. Reference Stevens, Zhong, Clercx, Ahlers and Lohse2009; Xia Reference Xia2011; Chong et al. Reference Chong, Wagner, Kaczorowski, Shishkina and Xia2018). Therefore, the LSC regime transition is expected to be accompanied by great alterations in heat transfer efficiency and turbulence intensity. The commensurate flow statistics are thus examined in the following for the fluid phase.
$(a)$
The mean temperature
$\langle T\rangle$
and
$(b)$
the Reynolds number
${\textit{Re}}$
versus
$\hat {\rho }_r$
at various
${\textit{Rr}}$
: (
),
${\textit{Rr}}=4.97\times 10^5$
; (
),
${\textit{Rr}}=1.49\times 10^6$
; (
),
${\textit{Rr}}=4.97\times 10^6$
; (
),
${\textit{Rr}}=1.49\times 10^7$
; (
),
${\textit{Rr}}=4.97\times 10^7$
; (
),
${\textit{Rr}}=1.49\times 10^8$
; (
),
${\textit{Rr}}=4.97\times 10^8$
. The red square, green triangle and blue circle superimposed on the lines represent, in sequence, the BFU, marginal and BFD regimes.
Similar to the IHC problem (Tasaka & Takeda Reference Tasaka and Takeda2005; Bouillaut et al. Reference Bouillaut, Lepot, Aumaître and Gallet2019; Wang et al. Reference Wang, Lohse and Shishkina2021), the stationary steady flow of interest is achieved when the total heat released by particles is balanced by that escaping from the top and bottom walls. Following Kuznetsov & Nield (Reference Kuznetsov and Nield2016), Wang et al. (Reference Wang, Lohse and Shishkina2021) and Creyssels (Reference Creyssels2021), the heat transfer efficiency and flow intensity are evaluated by the VT-averaged temperature
$\langle T \rangle$
and the Reynolds number defined as
$ Re=\sqrt {\langle \boldsymbol{u}^2\rangle } L/\nu$
, respectively. Note that
$1/\langle T \rangle$
takes the role of the Nusselt number and is regarded as an indicator of the heat transfer efficiency, because a lower
$\langle T \rangle$
is realized when the total heat released by particles is transferred more efficiently out of the PHC system via the necessary wall heat flux.
As indicated in figure 5, before the flow enters the BFD regime, both
$\langle T \rangle$
and
${\textit{Re}}$
decrease slightly with the increase of
$\hat {\rho }_r$
for all
${\textit{Rr}}$
. At a given
${\textit{Rr}}$
, lower
$\langle T \rangle$
corresponds to more efficient heat transfer for the PHC system (Tasaka & Takeda Reference Tasaka and Takeda2005; Kuznetsov & Nield Reference Kuznetsov and Nield2016; Bouillaut et al. Reference Bouillaut, Lepot, Aumaître and Gallet2019; Creyssels Reference Creyssels2021; Wang et al. Reference Wang, Lohse and Shishkina2021). For the BFU and marginal regimes (see figure 5
$a$
), the decrease of
$\langle T \rangle$
with increasing
${\textit{Rr}}$
indicates that more efficient heat transfer is facilitated by enhanced convective flow under stronger buoyancy. As expected, a corresponding increase in
${\textit{Re}}$
versus
${\textit{Rr}}$
is demonstrated in figure 5(
$b$
). When the flow regime transitions into the BFD regime (see figure 5
$a$
), the heat transfer efficiency demonstrates a complicated response to the increase of
$\hat {\rho }_r$
, namely, a great decrease of
$\langle T \rangle$
for
${\textit{Rr}}\lt 4.97\times 10^6$
while a slight decrease followed by a small increase for
${\textit{Rr}}\gt 4.97\times 10^7$
. Correspondingly,
${\textit{Re}}$
presents a substantial decrease versus
$\hat {\rho }_r$
, especially for
${\textit{Rr}}$
smaller than
$4.97\times 10^6$
. It indicates an overall suppression of the convective flow when the BFD regime is realized. The decrease of
$\langle T \rangle$
may be attributed to particle accumulation near the top and bottom walls, where the heat released by particles is expected to be transferred immediately and thus efficiently out of the system, before acting to heat the bulk fluid. This speculation will be further illustrated below by the HT-averaged temperature
$\overline {T}$
.
(
$a$
) The averaged temperature
$\langle T\rangle$
and
$(b)$
the Reynolds number
${\textit{Re}}$
versus
${\textit{Rr}}$
. The scaling laws of
$\langle T\rangle \sim {\textit{Rr}}^{-1/5}$
and
${\textit{Re}}\sim {\textit{Rr}}^{2/5}$
obtained from the IHC problem are plotted for comparison.
The above integral flow responses are closely related to the spatial distribution of heat-releasing particles. This can be examined here by the scaling behaviours of
$\langle T \rangle$
and
${\textit{Re}}$
via comparing with those of the IHC problem. As depicted in figure 6,
$\langle T \rangle$
and
${\textit{Re}}$
agree well with the scaling laws obtained for the IHC system (Goluskin & Van der Poel Reference Goluskin and Van der Poel2016; Goluskin Reference Goluskin2016; Wang et al. Reference Wang, Lohse and Shishkina2021), i.e.
$\langle T\rangle \sim {\textit{Rr}}^{-1/5}$
and
${\textit{Re}}\sim {\textit{Rr}}^{2/5}$
, respectively, at relatively small particle inertia of
$\hat {\rho }_r\leqslant 10$
. In other words, the heat transfer and flow intensity in the present PHC only slightly deviate from their IHC counterparts. Thus, it can be expected that the heat-releasing particles of small inertia may have a nearly uniform spatial distribution in the bulk flow. When the BFD flow regime is achieved, especially at high particle inertia (
$\hat {\rho }_r\geqslant 150$
),
$\langle T \rangle$
and
${\textit{Re}}$
deviate obviously from the
${-1/5}$
and
$2/5$
scaling laws. The spatial distribution of particles at high
$\hat {\rho }_r$
is thus implied to have a preferential concentration in the PHC, which is also a common situation in the particle-laden flows.
For the PHC problem, the wall heat transfer can also be measured by the thermal boundary layer thicknesses defined at the top
$(\lambda _t)$
and bottom
$(\lambda _b)$
walls (Wang et al. Reference Wang, Lohse and Shishkina2021; Du & Yang Reference Du and Yang2022),
\begin{equation} \lambda _{t, b}= \frac {\langle T \rangle }{|\textrm {d} \overline {T} /\textrm {d} {z} |_{{z}=1,0}}. \end{equation}
As seen in figure 7,
$\lambda _b$
is much greater than
$\lambda _t$
, indicating that relatively more heat is transferred out of the flow system through the top wall, and the thermal convection in the upper-half of the PHC is much stronger than that in the lower-half. This feature originates from the thermally unstable and stable stratification close to the top and bottom walls, respectively, similar to the IHC (Kuznetsov & Nield Reference Kuznetsov and Nield2016; Creyssels Reference Creyssels2021; Wang et al. Reference Wang, Lohse and Shishkina2021). It is basically different from the classic RBC, in which
$\lambda _b$
and
$\lambda _t$
are identical as a signature of intrinsic flow symmetry (Wondrak et al. Reference Wondrak, Sieger, Mitra, Schindler, Stefani, Vogt and Eckert2023; Wang et al. Reference Wang, Chen, Xu and Xi2024a
,
Reference Wang, Chen, Xu and Xib
). A comparison between figures 7(
$a$
) and 5(
$a$
) indicates that
$\lambda _t$
exhibits a similar
$\hat {\rho }_r$
dependence as
$\langle T\rangle$
for all
${\textit{Rr}}$
. Specifically,
$\lambda _t$
is slightly affected by increasing
$\hat {\rho }_r$
when the PHC remains in the BFU and marginal regimes. As the BFD regime is realized,
$\lambda _t$
decreases significantly for
${\textit{Rr}}\lt 4.97\times 10^7$
while changes subtly for
${\textit{Rr}}\gt 4.97\times 10^7$
versus
$\hat {\rho }_r$
. On the other hand,
$\lambda _b$
demonstrates a more intricate response to the change of
$\hat {\rho }_r$
when
${\textit{Rr}}$
is prescribed (see figure 7
$b$
). For the BFU and marginal flow regimes,
$\lambda _b$
decreases with the increase of
${\textit{Rr}}$
at a fixed
$\hat {\rho }_r$
, a similar response as that of
$\langle T \rangle$
discussed above. It undergoes a considerable decrease preceding a notable increase when the BFD regime is achieved at high
$\hat {\rho }_r$
, indicating that the heat transfer at the bottom wall is enhanced at first and then suppressed greatly during the flow regime transition.
Thermal boundary layer thicknesses close to
$(a)$
the top wall and
$(b)$
the bottom wall, denoted as
$\lambda _t$
and
$\lambda _b$
, respectively, as a function of
$\hat {\rho }_r$
at various
${\textit{Rr}}$
. Here, the lines and symbols are used following the exact styles of figure 5.
Representative flow quantities obtained for
${\textit{Rr}}=4.97\times 10^7$
at different
$\hat {\rho }_r$
:
$(a)$
the HT-averaged temperature
$\overline {T}$
and
$(b)$
the vertical convective heat flux
$\overline {{u_z} T}$
.
The local flow response to the influence of heat-releasing particles can be examined by the key variables obtained via HT averaging for the typical cases of
${\textit{Rr}}=4.97\times 10^7$
. As illustrated in figure 8(
$a$
), with the increase of
$\hat {\rho }_r$
, the flattened
$\overline {T}$
profiles in the bulk region demonstrate a non-monotonic variation manifested as a decrease (
$1\leqslant \hat {\rho }_r\leqslant 150$
) followed by a subsequent increase (
$150\leqslant \hat {\rho }_r\leqslant 1000$
). This behaviour is basically consistent with the integral flow response indicated by
$\langle T \rangle$
in figure 5(
$a$
). Notably, for the BFD regime especially at high
$\hat {\rho }_r$
, a stable thermal stratification quantified by the positive gradient of
$\overline {T}$
, appears in almost the entire fluid layer, except for the immediate vicinity of the top wall where an unstable one is maintained with a significant negative gradient of
$\overline {T}$
. These temperature profiles are expected to yield weak thermal convection near the bottom wall, while generating much more energetic ones in a very thin fluid layer near the top wall. This is confirmed by the thermal boundary layer thickness plotted in figure 7(
$a$
) and the HT-averaged vertical heat flux (
$\overline {{u_z} T}$
) in figure 8(
$b$
).
In fact,
$\overline {u_z T}$
exhibits a much larger positive and a small negative extremum in the top and bottom wall regions, respectively (see figure 8
$b$
). These two extremes appear due to the following flow behaviours: near the top wall, significant convective heat transfer is promoted by vigorous cold plumes that possibly penetrate down into the bulk region; near the bottom wall, convective flow is stirred up as warm plumes starting from the bulk and weakened in the stable thermal stratification. One should keep in mind that heat-releasing particles entrained by the thermal plumes or accumulated close to the walls may play an important role in the heat transfer efficiency. To this end, special interest is directed to the more energetic top wall region. For the BFD regimes, especially at higher
$\hat {\rho }_r$
(see figure 8
$a$
), a thinner fluid layer bedding over the
$\overline {T}$
maximum indicates preferential concentration of heat-releasing particles. Consequently, heat released by inertial particles is expected to escape immediately and thus more efficiently from the flow system through the top wall, before working to drive intense thermal convection in the cube. This expectation is confirmed by the remarkable decrease of
$\overline {u_z T}$
versus
$\hat {\rho }_r$
close to the top wall for the BFD regimes (see figure 8
$b$
). Moreover, the nearly linear profiles of
$\overline {T}$
and
$\overline {{u_z} T}$
in the bulk region together imply that heat sources (particle concentration) are almost vertically uniform there. This implication is inspired by our analytical findings about the balance between HT-averaged heat flux in non-uniform IHC (Chen et al. Reference Chen, Gao, Liao, Wan and Liu2024), and would be confirmed below by the analyses of particle concentrations.
Profiles of HT-averaged root-mean-square (r.m.s.) values of
$(a)$
temperature
$\overline {T}_{\textit{rms}}$
and
$(b)$
vertical velocity
$\overline {u}_{{z},rms}$
for
${\textit{Rr}}=4.97\times 10^7$
at different
$\hat {\rho }_r$
. Note that the vertical axis of panel (
$a$
) is scaled by
$10^{-3}$
.
The r.m.s. profiles of HT-averaged temperature and vertical velocity, i.e.
$\overline {T}_{\textit{rms}}$
and
$\overline {u}_{z,rms}$
, provide further details of the local flow features in response to the LSC regime transition. As viewed from figure 9(
$a$
),
$\overline {T}_{\textit{rms}}$
obtains higher peak values closer to the top wall with respect to the counterparts near the bottom wall for all
$\hat {\rho }_r$
. This observation indicates that more intense thermal fluctuations occur commensurate with higher heat flux across the top wall. It is consistent with the thinner thermal boundary layer of
$\lambda _t$
(see figure 7), as well as high peak values of
$\overline {T}$
and
$\overline {u_z T}$
(see figure 8
$a,b$
). During the flow regime transition, the peak values of
$\overline{T}_{\textit{rms}}$
exhibit a non-monotonic variation in the top and bottom wall regions, possibly due to the change in preferential concentration versus
$\hat {\rho }_r$
. As the BFD regime is realized, it is noted that
$\overline {T}_{\textit{rms}}$
presents a nearly linear profile in the bulk region (here,
$0.2\lt {z}\lt 0.8$
), which also has a weak but non-monotonic change with respect to
$\hat {\rho }_r$
.
As shown in figure 9(
$b$
), the maximum value of
$\overline {u}_{z,rms}$
occurs closer to the top wall in the BFU and marginal regimes (
$1\leqslant \hat {\rho }_r\leqslant 10$
), while shifting towards the bottom wall with substantial attenuation for the BFD regime at
$30\leqslant \hat {\rho }_r\leqslant 250$
. In particular, the profile of
$\overline {u}_{z,rms}$
at the highest
$\hat {\rho }_r=1000$
changes to resemble that of low
$\hat {\rho }_r$
belonging to the BFU and marginal regimes. Special interest is dedicated to the BFD regimes, which demonstrate a great enhancement of
$\overline {T}_{\textit{rms}}$
in the bulk region while demonstrating a considerable suppression in
$\overline {u}_{z,rms}$
(see figure 9
$a, b$
). The occurrence of preferential concentration of heat-releasing particles is further indicated in these BFD regimes, according to two previous illuminating findings. First, preferential concentration of heat-releasing particles near the top wall is found to hinder vertical convective flow (Tasaka & Takeda Reference Tasaka and Takeda2005). Second, heat-releasing particles that locally concentrate in the bulk region act to enhance the temperature fluctuation (Pouransari & Mani Reference Pouransari and Mani2018).
3.3. Particle distribution features
Heat source distribution has been known to greatly impact flow dynamics in thermal convection (Tasaka & Takeda Reference Tasaka and Takeda2005; Kuznetsov & Nield Reference Kuznetsov and Nield2016; Chen et al. Reference Chen, Gao, Liao, Wan and Liu2024). To this point, discussions on heat-releasing particle distribution are aimed to further interpret the above flow responses based on the typical cases of
${\textit{Rr}}=4.97\times 10^7$
. We attempt to substantiate the mechanistic understanding of the complex interplay between convective flows and particle dynamics behind the LSC regime transition. Following our recent work (Liao et al. Reference Liao, Chen, Wan, Liu and Lu2024), an Eulerian grid called ‘particle cells’ with a uniform spacing
$d_c$
larger than the particle diameter
$d_p$
is used to define the normalized particle concentration,
$c=N_{p,c} L^3/N_p d_c^3$
, with
$N_{p,c}$
being the number of particles positioned in a given particle cell. Here
$d_c=1/128$
corresponding to
$128^3$
particle cells in the cube is adopted to obtain converged results discussed below.
Particle distribution visualized by
$(a,b,c,d)$
time-averaged concentration and
$(e,f,g,h)$
instantaneous particle positions in typical flow regimes for
${\textit{Rr}}=4.97\times 10^7$
:
$(a,e)$
the BFU regime,
$(b,f)$
the marginal regime and
$(c,d)$
and
$(g,h)$
the BFD regimes.
The preferential concentration of heat-releasing particles is visualized by the time-averaged particle concentration
$\langle c \rangle _t$
and instantaneous particle positions in figure 10. For the BFU regime of
$\hat {\rho }_r=1$
, particles (heat sources) are observed to have a slight accumulation in the four top-wall corners, with a very small population indicated in the bottom-wall corners (see figure 10
$a,e$
). The remaining majority of heat-releasing particles have a nearly homogeneous distribution in the bulk and thus dominate the convective flow dynamics. It is consistent with empirical expectations that low-inertia particles generally have a homogeneous distribution. Consequently, the LSC is characterized as a BFU regime similar to the IHC situation (Goluskin & Spiegel Reference Goluskin and Spiegel2012). The local accumulation in these corners is possibly attributed to the interphase thermal coupling in the PHC. Namely, heat-releasing particles absent of gravity, once entering the top-wall corners, may be detained there by the weak circulation of surrounding warm plumes they generated.
For the marginal regime of
$\hat {\rho }_r=10$
, the preferential concentration is weakly enhanced in the corners intersected by the top/bottom and sidewalls (see figure 10
$b,f$
). As the BFD regime is realized, for example at
$\hat {\rho }_r=50$
, two distinct changes happen to the preferential concentration of particles. First, heat-releasing particles have a significant local accumulation, i.e. a heavy population under the entire top wall, and a small accumulation in the bottom corners (see figure 10
$c$
). Overall, this spatial change in particle distribution consistently accounts for the integral and local flow responses, as discussed in § 3.2. Second, the particle phase obtains structured clustering localized in the bulk and top/bottom-wall regions, as indicated by the instantaneous particle positions in figure 10(
$g,h$
). These two particle distribution features are confirmed and discussed in more detail in the following.
The vertical profiles of HT-averaged particle concentration
$\overline {c}$
for
${\textit{Rr}}=4.97\times 10^7$
.
Particle accumulation near the top and bottom walls can be quantified by the HT-averaged particle concentration
$\overline {c}$
for
${\textit{Rr}}=4.97\times 10^7$
, as depicted in figure 11. For the particle distribution over the bottom wall, an overturning change is demonstrated when the LSC transitions to the BFD regime. That is, the
$\overline {c}$
value increases over unit as
$\hat {\rho }_r$
increases, indicating that heat-releasing particles with lower inertia tend to migrate away from the bottom wall while high-inertia particles may accumulate there. Recalling the observations in figure 10, it is clear that higher
$\overline {c}$
near the top wall results from particles accumulating in the top corners for the BFU and marginal regimes (
$\hat {\rho }_r\leqslant 10$
). Meanwhile, particles have an approximately homogeneous distribution and consequently uniform
$\overline {c}$
almost throughout the bulk region. This confirms that linear profiles of
$\overline {T}$
and
$\overline {u_z T}$
indicate a relatively uniform distribution of heat sources (see figure 8). Of critical interest, all the BFD regimes are featured by significantly larger
$\overline {c}$
in the top-wall region, as compared with in the bottom-wall region. This feature becomes more evident with the increase of
$\hat {\rho }_r$
, accompanied by a gradual decrease of
$\overline {c}$
in the bulk. Thus, the profound preferential concentration occurring near the top wall is substantiated for the BFD regime. Clearly, particle accumulation near the top and bottom walls is responsible for the attenuation of mean temperature (
$\langle T \rangle$
) and convection intensity (
${\textit{Re}}$
). It is due to the fact that the heat released by particles there escapes immediately via conductive wall heat flux away from the cube, and therefore, does not work effectively to heat the bulk fluid and facilitate the thermal convection (Tasaka & Takeda Reference Tasaka and Takeda2005; Kuznetsov & Nield Reference Kuznetsov and Nield2016; Chen et al. Reference Chen, Gao, Liao, Wan and Liu2024).
Instantaneous particle positions (
$a{-}d$
) near the top wall (
${z}\gt 0.9$
) and (
$e{-}h$
) near the bottom wall (
${z}\lt 0.1$
) for
${\textit{Rr}}=4.97\times 10^7$
, topviewed typically for the BFU regime at
$\hat {\rho }_r=1$
, the marginal regime at
$\hat {\rho }_r=10$
and the BFD regimes at
$\hat {\rho }_r=50\,250$
.
The structured clustering of heat-releasing particles is visualized by the instantaneous particle positions in the top/bottom-wall regions at
${\textit{Rr}}=4.97\times 10^7$
(see figure 12). For the BFU (
$\hat {\rho }_r=1$
) and marginal (
$\hat {\rho }_r=10$
) regimes, the particle phase presents no evident structured clustering other than a randomly uniform distribution near the top and bottom walls. With the realization of the BFD regimes, structured clustering can be viewed clearly for the particles in the top/bottom-wall regions. Specifically, particle clustering is observed at the boundaries of well-organized cellular structures that are visible particularly near the top wall. This phenomenon is facilitated significantly by higher-inertia particles, for example, at
$\hat {\rho }_r=250$
(see figure 12
$d$
).
Interestingly, multiple cellular structures appear in the top-wall region, resembling the situation of the classic RBC where the flow domain is divided topologically into multiple convection cells. In the RBC problem, such multiple cellular structures are typically composed of energetic cold/hot plumes which can penetrate the entire thermal boundary layer (Wondrak et al. Reference Wondrak, Sieger, Mitra, Schindler, Stefani, Vogt and Eckert2023). It makes sense here that structured clustering of heat-releasing particles has a physical connection with the thermal plumes, which are distinctive in the top/bottom-wall regions for the BFD regimes. As evidence, cellular structures attain great differences close to the top and bottom wall regions (see figure 12
$c,g$
and
$d,h$
). This occurs due to distinctive convection behaviours near the top/bottom wall: near the top wall, vigorous cold plumes starting from the top wall intermittently penetrate into the bulk region, while near the bottom wall, weak warm plumes approach the bottom wall from the bulk. Such a mechanistic interpretation has been manifested by the flow statistics in § 3.2 and can be further acknowledged below.
Instantaneous particle positions and thermal plumes illustrated by arrows of particle velocity vector and contours of temperature fluctuations
$T^ \prime = T-\langle T \rangle _t$
, respectively, at
${\textit{Rr}}=4.97\times 10^7$
. The particles are extracted in the thin fluid layers immediately close to the top (
${z}\gt 0.99, a-d$
) and bottom (
${z}\lt 0.01, e-h$
) walls, while the temperature fluctuations are plotted on the planes (
${z}=0.99, a-d$
) and (
${z}=0.01, e-h$
), with
$(a,e)$
the BFU regime,
$(b,f)$
the marginal regime and
$(c,g)$
and
$(d,h)$
the BFD regimes.
The underlying mechanism of preferential concentration near the top/bottom walls is illustrated by instantaneous particle positions and thermal plumes in figure 13. Specifically, well-patterned cellular structures in the top-wall region are formed and identified by clustering lines of particle velocity vector arrows. However, they are absent in the bottom-wall region, especially for the BFU and marginal regimes. At low particle inertia (
$\hat {\rho }_r=1$
), these cellular structures have diverging particles at the centre and converging particles at the boundary, overlying well onto the convection cells which are centred at the rising warm (in red) plumes and bounded by the sinking cold (in blue) plumes (see figure 13
$a$
). Note here that such thermal convection cells are similar in topology to the classical Bénard cells (Gastine, Wicht & Aurnou Reference Gastine, Wicht and Aurnou2015). For the marginal regime of
$\hat {\rho }_r=10$
, a slight deviation of particle cellular structures from the convection cells of thermal plumes happens, indicating a lagging effect between the particle and fluid phases.
This lagging effect is enhanced in the BFD regime at
$\hat {\rho }_r=50$
in figure 13(
$c, g$
). It is manifested by a much more prominent deviation of particle cellular structures from thermal convection cells, i.e. notable deviation between the boundaries of these two flow patterns. For this BFD regime, an increase in number is indicated for the particle cellular structures and sinking cold plumes, along with a great suppression in strength for the rising warm plumes in the top-wall region (see figure 13
$c$
). The above aspects lead to an inefficient transport of particles away from the top wall, i.e. a near-top-wall preferential concentration. In contrast, both particle cellular structures and thermal convection cells are not well patterned in the bottom-wall region. It may be attributed to that the stable thermal stratification allows very few sinking warm plumes to survive and no rising cold plumes to be generated. Consequently, heat-releasing particles entering the bottom-wall region are trapped to have a relatively weaker aggregation (see also figure 11). It is consistent with the particle transport dynamics proposed by Du & Yang (Reference Du and Yang2022). See Supplementary Movies for more illustrations of particle transport dynamics.
As a further examination of particle aggregation, particle concentration obtained by VT averaging within a thin fluid layer close to the top (
$\langle c\rangle ^{t} =\langle c \rangle _{{z}\gt 0.9}$
) and bottom (
$\langle c\rangle ^{b} =\langle c \rangle _{{z}\lt 0.1}$
) walls, is manifested in figure 14. Overall, both
$\langle c\rangle ^{t}$
and
$\langle c\rangle ^{b}$
demonstrate a relatively small variation in the BFU regime. With the increase of
$\hat {\rho }_r$
,
$\langle c\rangle ^{t}$
obtains a monotonic increase for all
${\textit{Rr}}$
. Differently,
$\langle c\rangle ^{b}$
has an increase when the LSC regime changes from the BFU to the marginal regime, followed by a great decrease after reaching its peak value in the BFD regime. It agrees well with the typical prediction based on the HT-averaged particle concentration
$\overline {c}$
in figure 11, indicating that preferential concentration happens to all the BFD regimes. The variation of
$\langle c\rangle ^{t}$
(or
$\langle c\rangle ^{b}$
) presents a
$(\hat {\rho }_r,Rr)$
-dependent trend basically opposite to that of
$\lambda _t$
(or
$\lambda _b$
) depicted in figure 7. It confirms that the thermal boundary layer thickness is predominantly influenced by local accumulation of particles. Note that the occurrence of a peak value of
$\langle c\rangle ^{b}$
is also observed by Du & Yang (Reference Du and Yang2022) for particles of moderate inertia.
Particle concentration VT-averaged within a fluid layer of thickness 0.1 close to
$(a)$
the top and
$(b)$
the bottom walls, denoted as
$\langle c\rangle ^{t}$
and
$\langle c\rangle ^{b}$
, respectively, versus
$\hat {\rho }_r$
at different
${\textit{Rr}}$
. Here, the lines and symbols are used following the exact styles of figure 5.
To examine the intriguing impact of preferential concentration occurring near the sidewalls (see figure 10
$c,d$
), the sidewall particle concentration
$\widetilde {c}$
is examined in figure 15. It is indicated for
${\textit{Rr}}=4.97\times 10^7$
,
$\widetilde {c}$
has an identical value of approximately unit away from the sidewalls at
$\hat {\rho }_r=1, 3$
, confirming again that the BFU regime has a uniform particle distribution in the bulk region. The slight local accumulation approaching
$d=0$
may be caused by the centrifugal force analogous to the cases of high-inertia particles, although it is too weak to facilitate the LSC regime transition, as discussed below. When the flow transitions from the marginal to the BFD regime,
$\widetilde {c}$
demonstrates a gradual decline as
$d$
increases, followed by a notable rise (see figure 15
$b$
). It points to the fact that heat-releasing particles have an increasing accumulation in the sidewall regions, at the cost of decreasing population in the bulk. For
$\hat {\rho }_r\geqslant 150$
,
$\widetilde {c}$
changes to attain a decrease in the sidewall regions while a considerable increase to larger than unit in the bulk, indicating that preferential concentration is obtained in the regions of
$d\lt 0.03$
and
$d\gt 0.2$
. High-inertia particles are expelled by stronger centrifugal motion in the LSC flow, and prefer to accumulate near the non-slip and adiabatic sidewalls (
$d\lt 0.03$
) rather than in the bulk (
$d\gt 0.2$
), because they are decelerated when approaching the sidewalls and thus detained in the boundary layer. Therefore, the maximum value of
$\widetilde {c}$
is located near the sidewalls.
Sidewall particle concentration versus the minimum distance
$d$
to the sidewalls for
$ {\textit{Rr}}=4.97\times 10^7$
, illustrated as
$(a)$
linear plot and
$(b)$
logarithmic plot.
$(a)$
Mean temperature difference between the sidewalls and the centre
$(\delta T)$
,
$(b)$
mean vertical velocity of the particle phase near the sidewalls (
$\langle v_{{z}}\rangle ^s$
) versus
$\hat {\rho }_r$
at different
${\textit{Rr}}$
. Here, the lines and symbols are used following the exact styles of figure 5.
Considering the boundary conditions applied to the sidewalls, i.e. non-slip and adiabatic walls, particle (heat source) aggregation in the sidewall regions is deemed to raise local fluid temperature and in turn, the buoyancy force to facilitate the LSC regime transition. As a quantitative measure, the mean temperature difference between the sidewalls and the cube centre is introduced,
$\delta T= \langle T \rangle _{d=0}-\langle T \rangle _{d=0.5}$
, along with the mean vertical velocity of the particle phase near the sidewalls,
$\langle v_{{z}}\rangle ^s= \langle v_{{z}}\rangle _{d\lt 0.1}$
. It is indicated in figure 16
$(a)$
that the flow regime transition is characterized by a sign change of
$\delta T$
. Specifically,
$\delta T$
is negative in the BFU regimes and changes to be positive as the BFD regimes are achieved. The sign change of
$\delta T$
corresponds to a near-sidewall change from sinking cold plumes into rising warm plumes. As a consequence,
$\langle v_{{z}}\rangle ^s$
also undergoes a change in sign from negative to positive. Moreover,
$\langle v_{{z}}\rangle ^s$
experiences a decrease at the highest particle inertia
$\hat {\rho }_r=1000$
after peaking at a moderate
$\hat {\rho }_r$
for the BFD regime (see figure 16
$a$
). These findings demonstrate that the accumulation of heat-releasing particles near the sidewalls is responsible for the LSC regime transition.
3.4. Mechanical interpretation for LSC regime transition
A mechanical interpretation for the LSC regime transition is proposed as follows. In the BFU regime, particles of small inertia follow quite well the fluid motion and consequently have an almost uniform distribution in the cube, as evidenced in figure 11. When particle inertia (
$\hat {\rho }_r$
) is increased, the heat-releasing particles are prone to be expelled by the centrifugal force to aggregate near the sidewalls, giving birth to rising warm plumes that trigger an onset to the BFD flow regime when
$\hat {\rho }_r$
is large enough.
To this end, the centrifugal motion of particles within the LSC is pivotal in facilitating the transition from the BFU to the BFD regime, and quantitative characterization of particle transport is essential to verify the above mechanical interpretation and to identify the critical conditions.
Schematic diagram of transport dynamics of heat-releasing particles in the LSC structure of the BFU regime, with the red arrows denoting particle trajectory.
Following Eaton & Fessler (Reference Eaton and Fessler1994), the model described below is employed to analyse the transport mechanisms of inertial particles in the LSC structure of the BFU regime. As sketched in figure 17, the vortical flow in the LSC centre region is approximated by a conceptual Rankine vortex prescribed at the origin of polar coordinates (
$r,\theta$
). Given a constant angular velocity
$\varOmega$
(corresponding to a time scale
$1/\varOmega$
), the radial (
$u_r$
) and azimuthal (
$u_\theta$
) velocity components of the fluid phase are modelled as
\begin{equation} \begin{aligned} \begin{cases} u_r =0, \\ u_{\theta } =\varOmega r. \end{cases} \end{aligned} \end{equation}
In this context, the governing equations for a given particle at
$(r_p, \theta _p)$
are reduced as
\begin{equation} \begin{aligned} \begin{cases} \dfrac { {\textrm {d}} v_{r}}{ {\textrm {d}} t} = \dfrac {v_{\theta }^2}{r_p}+\dfrac {u_r-v_r}{\tau _p}, \\[12pt] \dfrac {{\textrm {d}} v_{\theta }}{{\textrm {d}} t} = \dfrac {u_{\theta }-v_{\theta }}{\tau _p}-\dfrac {u_{\theta }u_r}{r_p}.\\ \end{cases} \end{aligned} \end{equation}
The above formulae can be rewritten by
\begin{align*} \qquad\qquad\qquad\qquad \begin{cases} \ddot r_p +\dfrac {\dot r_p}{\tau _p} - r_p \dot \theta _p^2 &= 0, \qquad \qquad \quad\qquad\qquad\qquad\qquad \text{(3.5a)}\\[12pt] \ddot \theta _p + \left (\dfrac {1}{\tau _p} + 2\dfrac { \dot r_p}{r_p}\right ) \dot \theta _p &= \dfrac {\varOmega }{\tau _p}. \qquad \qquad\qquad\qquad\qquad\qquad \,\,\,\text{(3.5b)} \end{cases} \, \end{align*}
Clearly, particles released at any position (excluding those trapped right in the vortex centre) will eventually spiral away from the vortex core, with a positive radial velocity
$\dot r_p$
. According to (3.5b
), it is clear that the coefficient of the damping term,
$(({1}/{\tau _p}) + 2 ({ \dot r_p}/{r_p}))$
is always positive. This implies that the angular acceleration of particle
$\ddot \theta _p$
will approach zero as the system becomes fully developed. In this situation, the particle obtains a constant angular velocity
$\dot \theta _p$
, as well as a constant ratio between radial velocity and radial position,
${\dot r_p}/{r_p}$
.
In that, the differential equations (3.5) can be further reduced into algebraic equations by taking (
$\dot r_p /r_p, \dot \theta$
) as unknown variables, namely
\begin{equation} \begin{aligned} \begin{cases} \left (\dfrac {\dot r_p}{r_p} \right )^2 +\dfrac {1}{\tau _p}\dfrac {\dot r_p}{r_p}-\dot \theta _p^2 = 0,\\[12pt] \left (\dfrac {1}{\tau _p}+2\dfrac {\dot r_p}{r_p} \right )\dot \theta _p = \dfrac {\varOmega }{\tau _p}. \end{cases} \end{aligned} \end{equation}
For the above equations, the physically non-trivial solution reads as
\begin{equation} \begin{aligned} \begin{cases} \dfrac {\dot r_p}{r_p} =\varOmega \dfrac {f(\varOmega \tau _p) -1 +4 \varOmega ^2 \tau _p^2 f^3(\varOmega \tau _p)}{2\varOmega \tau _p },\\[10pt] \dot \theta _p = \varOmega f(\varOmega \tau _p), \end{cases} \end{aligned} \end{equation}
where
\begin{equation} f(x) = \left ( \dfrac {4\sqrt {x^2 + 1/16} -1}{8x^2} \right )^{1/2}. \end{equation}
Obviously, the particle radial velocity
$\dot r_p$
serves as a key indicator of the particle centrifugal motion driven by the vortex. Thus, an effective dimensionless parameter, defined as the ratio of the particle radial velocity to the vortical flow velocity,
$\varGamma =\dot r_p/(\varOmega r_p)$
, can be introduced to characterize the particle centrifugal motion in the LSC flows.
For the BFU regime of relatively small particle inertia (
$\varOmega \tau _p\ll 1$
), (3.7) yields
$\varGamma \sim \varOmega \tau _p$
. The angular velocity
$\varOmega$
is associated with the LSC structure by
$\varOmega \sim \sqrt {\langle \boldsymbol{u}^2\rangle }/L={\textit{Re}} \nu /L^2$
. With the definition of
$\tau _p$
given by (2.8),
$\varGamma$
can be estimated by
\begin{equation} \varGamma =\dfrac {\dot r_p}{\varOmega r_p} \sim \varOmega \tau _p \sim \dfrac {{\textit{Re}} \nu }{L^2} \dfrac {\rho _p d_p^2}{18\rho _f\nu } = {\textit{Re}} \hat {\rho }_r\dfrac {d^2_p}{18L^2}. \end{equation}
Recalling that the BFU regimes obtain a scaling fitting
${\textit{Re}} \sim {\textit{Rr}}^{2/5}$
in figure 6
$(b)$
, it leads to the following estimation:
\begin{equation} \varGamma \propto {\textit{Rr}}^{2/5} \hat {\rho }_r. \end{equation}
To quantify the mechanical interpretation proposed above,
$\varGamma$
obtained by (3.9) and (3.10) is plotted in figure 18 against
$\langle c \rangle ^s (=\langle c \rangle _{d\lt 0.1})$
, namely, the particle concentration averaged over particle cells near the sidewalls. Interestingly, the BFU and BFD regimes are well isolated in the phase space of (
${\textit{Re}}\hat {\rho }_r d_p^2/18L^2$
,
$\langle c \rangle ^s$
) by a small zone of the marginal regime (see figure 18
$a$
). A similar observation is obtained in figure 18(
$b$
) for the phase space of (
${\textit{Rr}}^{2/5} \hat {\rho }_r$
,
$\langle c \rangle ^s$
), confirming that the particle centrifugal motion plays a driving role for the LSC regime transition. In this sense, the onset condition of the LSC regime transition can be roughly estimated by (
${\textit{Rr}}^{2/5} \hat {\rho }_r$
,
$\langle c \rangle ^s$
) values. Specifically, when
${\textit{Rr}}^{2/5} \hat {\rho }_r$
is larger than
$5\times 10^4$
or
$\langle c \rangle ^s$
exceeds approximately
$3.5$
, the LSC structure is likely to transition from the BFU regime to the BFD regime. Note that the marginal regime corresponds to a narrow range of
${\textit{Rr}}^{2/5} \hat {\rho }_r$
. It thus leads to a reasonable explanation of the scaling of
$\hat {\rho }_r \propto {\textit{Rr}}^{-2/5}$
fitted for the marginal regime in figure 4.
We are aware that the above analysis starts from the conceptual vortex–particle model, and the LSC structure is modelled simply without taking into account the constraint of solid walls. Nevertheless, the key features of centrifugal-driven particle transport in the LSC flow are identified as the mechanistic interpretations.
4. Conclusion
Thermal convection driven by heat-releasing particles in a cube is studied via four-way coupling DNS based on the DEM approach, examining the influence of the thermal buoyancy (
${\textit{Rr}}$
) and particles inertia (
$\hat {\rho }_r$
) on the LSC structure dynamics. Wide parameter ranges of
${\textit{Rr}}=4.97\times 10^{5} - 4.97 \times 10^{8}$
and
$\hat {\rho }_r=1- 1000$
are considered along with
$Pr=1$
, to interpret the underlying correlations between the LSC regime transition and the commensurate flow statistics, and the particle distribution features.
An intriguing flow transition has been found via the identification of distinct LSC structures, involving in sequence three typical LSC regimes, i.e. the BFU regime, the marginal regime and the BFD regime. The BFU regime emerging typically at lower values of
$\hat {\rho }_r$
, has rising flow in the bulk and sinking flow close to the sidewalls. It is manifested as two dominating, symmetrical and counter-rotating roll cells in the vertical diagonal plane. As
$\hat {\rho }_r$
is increased to a moderate value, the flow transitions into the marginal regime which is characterized by a bulk composed of both rising-up and sinking-down fluids, and the LSC is manifested as a very large, highly distorted roll cell in the vertical diagonal plane. Further increasing
$\hat {\rho }_r$
leads to the onset of the BFD regime, which obtains a flow topology similar to the BFU regime but with an opposite LSC flow. Two first-order central moments of time-averaged vertical velocity (
$M_x, {M_y}$
) are employed as the diagnostic indicators to quantify the flow regimes. It is revealed that particle inertia plays a critical role in triggering the LSC regime transition, and the marginal regime exhibits a power-law scaling of
$-2/5$
in the (
${\textit{Rr}},\hat {\rho }_r$
) phase diagram.
The influence of the LSC regime transition on the heat and momentum transfer is elucidated by key flow statistics of the fluid phase. The (
${\textit{Rr}},\hat {\rho }_r$
) dependence of the mean temperature
$\langle T \rangle$
and Reynolds number
${\textit{Re}}$
indicates a comprehensive integral response of the heat transfer efficiency and flow intensity. It is found that when particle inertia is relatively small,
$\langle T \rangle$
and
${\textit{Re}}$
follow the scaling laws analogous to those observed in the IHC problem, specifically
$\langle T \rangle \sim {\textit{Rr}}^{-1/5}$
and
${\textit{Re}} \sim {\textit{Rr}}^{2/5}$
. A significant decline in
$\langle T \rangle$
and
${\textit{Re}}$
takes place as
$\hat {\rho }_r$
increases and the flow regime transitions into the BFD regime. The thermal boundary layer thicknesses at the top and bottom walls,
$\lambda _t$
and
$\lambda _b$
, exhibits similar responses to those of
$\langle T \rangle$
to the change of
${\textit{Rr}}$
. The HT-averaged profiles of temperature
$\overline {T}$
, vertical convective heat flux
$\overline {{u_z} T}$
, temperature fluctuations
$\overline {T}_{\textit{rms}}$
and velocity fluctuations
$\overline {u}_{{z},rms}$
point to an overall suppression of thermal convection when the flow enters into the BFD regime.
Particle distribution is analysed based on typical flow regimes to substantiate the mechanistic understanding of the complex interplay between the LSC flow and particle dynamics. Significant local accumulation occurs as
$\hat {\rho }_r$
increases, leading to great modification in the flow dynamics. Specifically, particles aggregating near the sidewalls heat the local surrounding fluid, inducing rising warmer plumes that in turn account for the LSC regime transition. On the other hand, particles accumulating greatly near the top/bottom walls result in attenuated thermal convection in the BFD regime. Intriguingly, the instantaneous particle distribution demonstrates well-patterned cellular structures in the top-wall region. These cellular structures overlay well on the thermal convection cells for small
$\hat {\rho }_r$
, but become notably deviated for the BFD regimes due to the lagging effect enhanced by higher
$\hat {\rho }_r$
. These observations demonstrate the physical connection between the structured clustering of heat-releasing particles and the thermal plumes especially in the top wall region.
A mechanistic interpretation for the LSC regime transition is substantiated based on a conceptual vortex–particle model. Specifically, a simplified analytical description is obtained for the particle transport dynamics, and an effective dimensionless parameter
$\varGamma$
is proposed to characterize the particle centrifugal motion in the LSC flow. This presents an easy estimation as a scaling law
$\varGamma \sim {\textit{Rr}}^{2/5}\hat {\rho }_r$
, which confirms along with the particle concentration near the sidewalls that that the centrifugal motion of heat-releasing particles plays a driving role in the LSC regime transition.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.10577.
Acknowledgements
The authors thank Dr Y.-H. Du for his instructive discussions. The simulations were performed on the Supercomputing Center of the University of Science and Technology of China and the Hefei Advanced Computing Center.
Funding
This work is supported by the Science Challenge Project; the National Natural Science Foundation of China under grant nos. 12388101, 92252202, 12172353, 12402257 and 12322210; the Fundamental Research Funds for the Central Universities and the USTC Research Funds of the Double First-Class Initiative.
Declaration of interests
The authors report no conflict of interest.
Appendix A. More details for DEM simulation
Following previous work (Schmeeckle Reference Schmeeckle2014; Zhu et al. Reference Zhu, Zhou, Yang and Yu2007b
), particle–particle interactions are formulated by the extensively applied soft-sphere approach taking advantage of its high accuracy. The collision force
$\boldsymbol{F}_{\textit{ij}}^{c}$
between a pair of contacting particles
$i$
and
$j$
can be evaluated by using an elastic spring and a viscous dashpot, having its normal (
$\boldsymbol{F}_{\textit{ij,n}}$
) and tangential (
$\boldsymbol{F}_{\textit{ij,t}}$
) components calculated as
\begin{align} \boldsymbol{F}_{\textit{ij}}^{c}&=\boldsymbol{F}_{\textit{ij,n}}+\boldsymbol{F}_{\textit{ij,t}}, \end{align}
\begin{align} \boldsymbol{F}_{\textit{ij,n}} &= k_n{\boldsymbol \delta }_n - \eta _n \boldsymbol{v}_n, \end{align}
\begin{align} \boldsymbol{F}_{\textit{ij,t}} &= \begin{cases} K_t \delta _t -\eta _t \boldsymbol{v}_t,& \,\boldsymbol{F}_{\textit{ij,t}}\lt \unicode{x03BC} _c \boldsymbol{F}_{\textit{ij,n}}, \\ \unicode{x03BC} _c \boldsymbol{F}_{\textit{ij,n}}, & \,\boldsymbol{F}_{\textit{ij,t}} \geqslant \unicode{x03BC} _c \boldsymbol{F}_{\textit{ij,n}}. \end{cases} \end{align}
For the contacting particles
$i$
and
$j$
,
$\delta _n$
and
$ \delta _t$
are the normal and tangential components of relative displacement, respectively;
$\boldsymbol{v}_n$
and
$\boldsymbol{v}_t$
are the normal and tangential components of relative velocity, respectively;
$k_n$
and
$k_t$
are the normal and tangential elastic coefficients, respectively;
$\eta _n$
and
$\eta _t$
the normal and tangential viscoelastic damping coefficients, respectively. Here
$\unicode{x03BC} _c$
is the friction coefficient between contacting particles.
The above coefficients of the soft-sphere approach are calculated by the Hertz contact model involving the material properties, given as
\begin{align} k_n &= \frac {2Y}{3(1-\nu _p^2)}\sqrt {\delta _n d_p/4}, \end{align}
\begin{align} k_t &= \frac {2Y}{(2-\nu _p)(1+\nu )}\sqrt {\delta _n d_p/4}, \end{align}
\begin{align} \eta _n &= -2\sqrt {\frac {5}{6}} \frac {\ln {e}}{\sqrt {\ln ^2{e}+\pi ^2}} \sqrt {S_n m}, \end{align}
\begin{align} \eta _t &= -2\sqrt {\frac {5}{6}} \frac {\ln {e}}{\sqrt {\ln ^2{e}+\pi ^2}} \sqrt {S_t m}, \end{align}
\begin{align} S_n &= \frac {Y}{1-\nu _p^2}\sqrt {\delta _n d_p/4}, \end{align}
\begin{align} S_t &= \frac {2Y}{(2-\nu _p)(1+\nu _p)}\sqrt {\delta _n d_p/4}, \end{align}
where
$Y$
,
$G$
,
$\nu _p$
and
$e$
represent in sequence the Young’s modulus, shear modulus, Poisson ratio and coefficient of restitution of the particle material;
$\delta _n$
is the magnitude of the normal displacement
$\delta _n$
between two contacting particles.
Dimensional parameter settings.
Validations based on the mean temperature
$(a)$
and vertical velocity
$(b)$
at
$y=0.5$
and
$z=0.5$
, compared with Gereltbyamba & Lee (Reference Gereltbyamba and Lee2019) for
$d_p=75\,{\rm \unicode{x03BC} m}$
and
$\varPhi _v=10^{-4}$
. Here,
$x$
denotes the wall-normal direction of the flow domain.
Verifications based on the HT-averaged temperature
$(a)$
and particle concentration
$(b)$
at
${\textit{Rr}}=4.97 \times 10^8$
.
The conductive heat
$Q_{\textit{ij}}^{c}$
exerted on particle
$i$
by particle
$j$
are modelled as (Li & Mason Reference Li and Mason2000; Bodhisattwa et al. Reference Bodhisattwa, Fernando and Silvina2006)
\begin{align} Q_{\textit{ij}}^{c} &= h_{c,ij}(T_j-T_i) , \end{align}
\begin{align} h_{c,ij} &= 2\kappa _p \sqrt {A_c} , \end{align}
where
$h_{c,\,ij}$
is the heat transfer coefficient between particles
$i$
and
$j$
, depending on the particle thermal conductivity
$\kappa _p$
and the particle contact area
$A_c$
. Here
$A_c$
is calculated by the particle diameter
$d_p$
and the normal displacement
$\delta _n$
between two contacting particles. More details are given in Brilliantov et al. (Reference Brilliantov, Spahn, Hertzsch and Pöschel1996) and Ai et al. (Reference Ai, Chen, Rotter and Ooi2011).
Neglecting the short-range perturbations around the particle, the momentum and thermal source terms, i.e.
$ \boldsymbol{F}_{\kern-1.5pt f}$
and
$Q_{\kern-1.5pt f}$
in (2.1), are modelled by the Dirac distributions
$\delta (\boldsymbol{x})$
,
\begin{align} \boldsymbol{F}_{\kern-1.5pt f} &= - \sum _{j}^{N_p} \frac {C_d {\textit{Re}}_p}{St_p}(\boldsymbol{u}_i -\boldsymbol{v}_i) \delta (\boldsymbol{x}-\boldsymbol{x}_j), \end{align}
\begin{align} Q_{\kern-1.5pt f} &= - \sum _{j}^{N_p} \frac {{\textit{Nu}}_p}{St_{\theta }}(T_{\textit{fi}} -T_i) \delta (\boldsymbol{x}-\boldsymbol{x}_j), \end{align}
where
$\boldsymbol{x}_j$
denotes the position of particle
$j$
. More details about the deposition of
$ \boldsymbol{F}_{\kern-1.5pt f}$
and
$Q_{\kern-1.5pt f}$
are described in Zhu et al. (Reference Zhu, Zhou, Yang and Yu2007b
) and Li & Mason (Reference Li and Mason2000).
All the dimensionless parameters involved are defined based on the parameter settings listed below in table 1 for the fluid flow, as well as the material properties of the particles. For the purpose of scientific interest, the parameters we set here are not strictly tailored to a real problem but from a dimensionless perspective.
Appendix B. Validation and verification
Although the present DNS-DEM approach has been successfully used for a wide range of applications (Xie et al. Reference Xie, Hu, Zhu, Yu and Pähtz2023a , Reference Xie, Zhu, Hu, Yu and Panb ; Schmeeckle Reference Schmeeckle2014), extensive simulations are conducted for further validation. To this end, the particle-laden flows in a differentially heated cubic cavity are solved to compare with the benchmark work of Gereltbyamba & Lee (Reference Gereltbyamba and Lee2019), using the two-way coupling model for the momentum and thermal exchange between the particle and fluid phases. As indicated in figure 19, good agreements are achieved for the typical results including the time-averaged temperature and vertical velocity. They assure that the present DNS-DEM solver works reliably to produce the essential flow features of the particle-laden thermal convection.
We have demonstrated that the grid resolution used is adequate to resolve the smallest turbulence scales
$(\eta _u,\, \eta _T)$
. As a further verification, test simulations are performed using two different meshes:
$N_x\times N_y\times N_z=120 \times { 120 \times 160}$
(Mesh1) and
$160\times {160 \times 384}$
(Mesh2), with
${\textit{Rr}}=4.97\times 10^8$
and
$\hat {\rho }_r=1$
,
$1000$
. In figure 20, the averaged temperature profiles and particle concentrations demonstrate negligible differences for the Mesh1 and Mesh2 cases, thus confirming the grid convergence.
Appendix C. Numerical details
Numerical details and key response parameters of each case are presented in table 2.
Summary of computational details and key response parameters.
