2.1 Working assumptions and dimensional equations

We consider the equation for the evolution of the fluid velocity field $\boldsymbol{u}$ coupled to the temperature field $T$ through the Oberbeck–Boussinesq approximation in a tri-dimensional periodic domain of linear size $H$ . The fluid temperature variations are proportional to density variations via the thermal expansion coefficient $\unicode[STIX]{x1D6FC}$ : $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}_{0}(1-\unicode[STIX]{x1D6FC}(T-T_{0}))$ in the low Mach number limit, where $\unicode[STIX]{x1D70C}_{0}$ and $T_{0}$ are the reference fluid density and temperature. The continuity equation and the momentum equation are:

with $\boldsymbol{u}$ the fluid velocity, $D_{t}=\unicode[STIX]{x2202}_{t}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ the material derivative, $\unicode[STIX]{x1D708}$ the kinematic viscosity, $p$ the pressure field, $g$ the gravity (in the $z$ direction) and $\boldsymbol{f}=\boldsymbol{F}/\unicode[STIX]{x1D70C}_{0}$ accounts for the exchange of momentum with the disperse phase. The temperature equation is:

with $\unicode[STIX]{x1D705}$ the thermal diffusivity and $q=Q/\unicode[STIX]{x1D70C}_{0}c_{f}$ , where $Q$ is the thermal source term per unit volume due to the heat exchanged with the disperse phase and $c_{f}$ is the fluid heat capacity.

The evolution equations for the particle velocity, position and temperature are written in the Lagrangian framework. We consider spherical particles much denser than the fluid and much smaller than the minimum scale of the flow. Retaining only the inertia, the Stokes drag and the gravitational force, the particle equations of motion, in the point particle approximation, are:

where $\boldsymbol{x}_{p}$ is the particle position coordinate, $\boldsymbol{u}_{p}$ is the particle velocity, $\unicode[STIX]{x1D70C}_{p}$ is the particle density and $\unicode[STIX]{x1D70F}_{p}=\unicode[STIX]{x1D70C}_{p}d_{p}^{2}/\unicode[STIX]{x1D70C}_{0}18\unicode[STIX]{x1D708}$ is the particle relaxation time, with $d_{p}$ the particle diameter. Considering the heat exchanged with the fluid and an external heat source, the temperature equation for a particle reads:

where $T_{p}$ is the particle temperature, $m_{p}$ is the mass of a particle, $c_{p}$ is the particle heat capacity and $\unicode[STIX]{x1D70F}_{th}=(3/2)(\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705})(c_{p}/c_{f})\unicode[STIX]{x1D70F}_{p}$ is the thermal relaxation time. $\unicode[STIX]{x1D6F7}_{p}$ is the external heat flux received by one particle and will be assumed constant throughout the paper. We consider small particle loading, such that particle collisions can be ignored. Thermophoresis, i.e. the transport of suspended particles in response to gradients in the fluid temperature, might also, in principle, play a role in the considered physical scenario. However, thermophoresis is only expected to be significant for particles of extremely low inertia, and outside of the range of parameters considered here. Its investigation is beyond the scope of the present study.

Neglecting the short-range perturbations around the particle, the momentum and thermal source terms $F$ and $Q$ in (2.2) and (2.3) are formally taken into account by Dirac distributions of the velocity and temperature differences:

In order to have a zero-mean forcing term in (2.2) and (2.3) the reference temperature is defined as the mean fluid temperature $T_{0}=\overline{T}$ where the $\overline{\bullet }$ denotes spatial averaging. From (2.3), accounting for the homogeneity of the flow, the mean temperature rate of change $\unicode[STIX]{x1D6FD}$ is expressed as: $\unicode[STIX]{x1D6FD}=d_{t}\overline{T}=\overline{q}$ . Introducing the temperature fluctuation around the spatially averaged fluid temperature $\unicode[STIX]{x1D703}=T-\overline{T}$ and $\unicode[STIX]{x1D703}_{p}=T_{p}-\overline{T}$ , (2.2), (2.3) and (2.11) become:

with $q^{\prime }=q-\unicode[STIX]{x1D6FD}$ . We will consider that the system is in a quasi-stationary state assuming that the characteristic time scale of the heating, $T_{0}/\unicode[STIX]{x1D6FD}$ , is large compared to the relevant time scale of the dynamics.

Using (2.6) and (2.8), one obtains the following relation between $\unicode[STIX]{x1D6FD}$ and $\unicode[STIX]{x1D6F7}_{p}$ :

where $\unicode[STIX]{x1D712}=\unicode[STIX]{x1D6FC}_{v}\unicode[STIX]{x1D70C}_{p}c_{p}/\unicode[STIX]{x1D70C}_{0}c_{f}$ with $\unicode[STIX]{x1D6FC}_{v}=\overline{n}m_{p}/\unicode[STIX]{x1D70C}_{p}$ the volume fraction of the disperse phase and $\overline{n}=N_{p}/H^{3}$ the particle number density. In the second term on the right-hand side of (2.12) we introduce the average over the set of particles $\langle \bullet \rangle _{p}$ . This term represents the increase of the thermal energy of the disperse phase relatively to the fluid phase. After a transient time of the order of $\unicode[STIX]{x1D70F}_{th}$ , this term vanishes in the statistically steady-state system.

In appendix A the equations (2.1), (2.9) and (2.10) are derived from the compressible Navier–Stokes equations assuming a small enough heat source term along with the usual assumptions required in the Oberbeck–Boussinesq approximation to be valid (Spiegel & Veronis Reference Spiegel and Veronis1960; Génieys & Massot Reference Génieys and Massot2001; Dumont et al. Reference Dumont, Genieys, Massot and Volpert2002; Shirgaonkar & Lele Reference Shirgaonkar and Lele2006), namely a vanishingly small Mach number and a small domain size in comparison to the length scale of the hydrostatic pressure variation ( $H\ll L_{\unicode[STIX]{x1D70C}}=P_{0}/\unicode[STIX]{x1D70C}_{0}g$ , $P_{0}$ being a reference pressure).

As an example application Zamansky et al. (Reference Zamansky, Coletti, Massot and Mani2014) considered the case of radiative heating of a particle–fluid mixture. Considering a transparent fluid and assuming an optically thin medium, the heat received by each particle is $\unicode[STIX]{x1D6F7}_{p}=\unicode[STIX]{x03C0}/4d_{p}^{2}\unicode[STIX]{x1D6F7}$ with $\unicode[STIX]{x1D6F7}$ the radiative heat flux density.

2.2 Fast thermal relaxation

We will mostly focus on the case of negligibly small particle thermal inertia ( $c_{p}/\!c_{f}\,\rightarrow \,0$ or equivalently $\unicode[STIX]{x1D70F}_{th}\rightarrow 0$ ). In this limit, particles are in thermal equilibrium with the surrounding fluid and (2.11) is replaced by:

With this assumption, the heat exchanged between the particles and the fluid is equal to the energy flux received by the particles, and (2.12) simplifies to:

and (2.8) becomes:

and one can express $q^{\prime }$ as:

2.3 Kinetic energy and thermal fluctuation budgets

From (2.9) the rate of change of the turbulent kinetic energy $K=\boldsymbol{u}^{2}/2$ is expressed as:

where $\unicode[STIX]{x1D700}=\unicode[STIX]{x1D708}(\unicode[STIX]{x1D735}u)^{2}=\unicode[STIX]{x1D708}\unicode[STIX]{x2202}_{x_{j}}u_{i}\unicode[STIX]{x2202}_{x_{i}}u_{j}$ is the rate of turbulent kinetic energy dissipation and the last two terms on the right-hand side represent the production of $K$ due to buoyancy (and the rate of variation of the potential energy), with $w$ the vertical component of the velocity and the variation of $K$ caused by the coupling with the dispersed phase. Assuming stationarity and neglecting the momentum coupling between phases gives:

From (2.10) the temperature variance budget gives:

where $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x1D705}(\unicode[STIX]{x1D735}\unicode[STIX]{x1D703})^{2}$ is the rate of dissipation of the temperature variance. The production term $\unicode[STIX]{x1D703}q^{\prime }$ can be re-expressed with (2.16), in the case of vanishing particle thermal inertia, as:

The fluid temperature at the particle position $\unicode[STIX]{x1D703}_{p}$ is, in this case:

with $V$ the volume of the fluid domain. Then assuming the statistical stationarity in (2.19), one obtains:

where $\overline{\unicode[STIX]{x1D703}}=0$ has been used.

2.4 Characteristic scales and non-dimensional equations

In the present setting, the relevant scales of the dynamics are not imposed by the forcing, but result from the inter-phase coupling. We estimate characteristic scales of the fluid motion in response to the local heating as follows. We assume that at those scales, conductive and convective heat transfer are balanced, i.e. a unit Rayleigh number:

and that inertia and viscous forces are also balanced i.e. a unit Reynolds number:

where $\unicode[STIX]{x1D703}_{\ast }$ , $\ell _{\ast }$ , $u_{\ast }=\ell _{\ast }/t_{\ast }$ and $t_{\ast }$ are the characteristic temperature, length, velocity and time of the flow at small scales. The characteristic heating rate of the system $\unicode[STIX]{x1D6FD}$ relates the temperature scale to the time scale:

Equations (2.23)–(2.25) give for the characteristic temporal and length scales:

As will be shown in § 4.3, the results of the simulations indicate that these scales are indeed of the order of the dissipative scales of the flow.

The following set of non-dimensional parameters is then introduced:

1. (a) the Stokes number: $St=\unicode[STIX]{x1D70F}_{p}/t_{\ast }$ ,

2. (b) the non-dimensional particle number density: $C=\overline{n}\ell _{\ast }^{3}$ (i.e. the average number of particles in a volume $\ell _{\ast }^{3}$ ),

3. (c) the non-dimensional domain size: $\unicode[STIX]{x1D6FE}=H/\ell _{\ast }$ ,

4. (d) the Froude number: $Fr=(g(1-\unicode[STIX]{x1D70C}_{0}/\unicode[STIX]{x1D70C}_{p})t_{\ast }^{2}/\ell _{\ast })^{-1/2}$ (i.e. the ratio of the particle gravitational acceleration and the buoyancy-induced fluid acceleration),

5. (e) the Prandtl number: $Pr=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$ ,

6. (f) the density ratio: $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{0}$ ,

7. (g) the specific heat capacity ratio $c_{p}/c_{f}$ .

Note that the non-dimensional particle number density is chosen over the particle volume fraction because the flow is influenced more by the distribution of the heat sources than by their volume since we consider particles with a vanishingly small diameter. However the following expression gives the volume fraction $\unicode[STIX]{x1D6FC}_{v}$ in terms of the other non-dimensional parameters:

The (2.1), (2.4), (2.5), (2.9) (2.10) and (2.16) non-dimensionalized by $t_{\ast }$ , $\ell _{\ast }$ , $u_{\ast }$ and $\unicode[STIX]{x1D703}_{\ast }$ can be re-written, in the fast thermal relaxation limit, and neglecting the reactions of the particles on the fluid momentum:

The thermal source term in (2.31) is:

where $\unicode[STIX]{x1D6FF}$ denotes the non-dimensional Dirac distribution.

As seen from (2.30) the absence of non-dimensional parameters in the momentum equation reflects the balance between inertia, buoyancy and viscous forces on which are based the scales $\ell _{\ast }$ and $t_{\ast }$ . On the other hand, the particles dynamics is characterized by two non-dimensional parameters, $St$ and $Fr$ , which control the overall dynamics of the system together with the particle number density. Finally $\unicode[STIX]{x1D6FE}$ reflects the boundary conditions delimiting the domain size.

3.1 Parameter values

Tables 1 and 2 give the parameters of the different simulations performed. In a first set of simulations, we consider the effect of varying the Stokes number of the particles neglecting their thermal inertia, direct feedback on the fluid momentum and gravitational settling. As seen in table 1, different values of $\unicode[STIX]{x1D6FE}$ (at constant $C$ ) and different values of $C$ (at constant $\unicode[STIX]{x1D6FE}$ ) are considered. The sets of simulations presented in table 2 were carried out to analyse the influence of the particle thermal inertia, momentum two-way coupling, as well as the effect of the particle settling. Note that for all simulations presented in tables 1 and 2, we set $Pr=1$ . As seen in the tables, the mesh resolution $\unicode[STIX]{x1D6E5}/\ell _{\ast }=\unicode[STIX]{x1D6FE}/N$ ranges from 0.55 to 0.75 which is enough to resolve the dissipative scale.

As in Zamansky et al. (Reference Zamansky, Coletti, Massot and Mani2014), all the simulations are initialized with randomly distributed particles and with zero velocity and temperature fluctuations for both fluid and particles. After a transient state, the flow reaches a statistically steady state which is our interest in this paper.

Table 1.

Parameters of the different simulations with zero particle thermal inertia and settling velocity. $N$ is the size of the mesh in each direction, $N_{p}$ is the number of particles. $St=0.003\sim 29.36$ refers to $St=0.003$ , 0.019, 0.074, 0.352, 1.064, 7.343 and 29.36, and MTWC indicates if the momentum two-way coupling $f$ is taken into account in (2.2).

Table 2.

Parameters of the different simulations for finite particle thermal inertia and settling velocity. For details see the caption of table 1. Here $0.003\sim 29.36^{a}$ refers to $St=0.003$ , 0.019, 0.074, 0.352, 1.064, 4.26, 7.343, 12.69 and 29.36. Additionally for all simulations $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{0}=909$ .

4.1 Clustering transition

In figure 1 we present visualizations of the instantaneous particle concentration field and fluid temperature fluctuations $\unicode[STIX]{x1D703}$ in a vertical plane for three levels of particle inertia: $St=0.003$ , $St=0.352$ and $St=7.343$ . In this figure the concentration is normalized by the average particle concentration $\overline{n}$ and $\unicode[STIX]{x1D703}$ is normalized by $\overline{\unicode[STIX]{x1D703}^{2}}^{1/2}$ . It is observed that for the smallest particle inertia, both fields are fairly homogeneous with only small-scale patterns in the temperature field. At the intermediate Stokes number, highly concentrated clusters of particles appear along with a more coherent distribution the temperature field organized in sheet-like structures. It can be visually checked that, for this Stokes number, the largest temperature fluctuations are highly correlated with the high particle concentration spots. This is in agreement with the conclusions drawn previously in Zamansky et al. (Reference Zamansky, Coletti, Massot and Mani2014). For the high Stokes number, the particles distribution becomes again much more homogeneous although very large areas, commensurate with the size of the overall domain, depleted of particles subsist. The temperature field also present large regions of positive and negative fluctuations with typical mixing layer structures in between.

Figure 1.

Visualization in a $x{-}z$ plane: (a,c,e) of the logarithm of the particle concentration field (as computed by (3.1)) normalized by the mean concentration $n/\overline{n}$ , and (b,d,f) the fluid temperature fluctuations normalized by its standard deviation $\unicode[STIX]{x1D703}/\overline{\unicode[STIX]{x1D703}^{2}}^{1/2}$ . For $\unicode[STIX]{x1D6FE}=570$ and $C=0.19$ , at three Stokes numbers: $St=0.003$ (a,b), $St=0.352$ (c,d) and $St=7.343$ (e,f).

Figure 2 depicts the structure of the turbulent kinetic energy dissipation rate $\unicode[STIX]{x1D700}$ and of the scalar dissipation rate $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}$ , for the same snapshots as in figure 1. Both quantities are normalized by their respective mean values. The effect of the Stokes number on the spatial organization of the turbulence dissipation rate is not very pronounced. Filaments of intense dissipation are observed for the three particle inertia, although for the two largest Stokes numbers the filaments appear to be more intense and corrugated. The temperature dissipation rate is fairly homogeneous for the smallest particle inertia. No structures are seen and only small deviations from the mean value are observed. At $St=0.352$ and $St=7.343$ instead, the temperature dissipation rate presents filaments of very high intensity (more than two decades larger than the mean values, as will be confirmed in figure 18). By comparing with the temperature fields in figure 1 we observe, as expected, that the filaments are found where the temperature gradient is large.

Figure 2.

Visualization in a $x{-}z$ plane: (a,c,e) logarithm of the turbulent kinetic energy dissipation rate normalized by its mean value $\log _{10}(\unicode[STIX]{x1D700}/\overline{\unicode[STIX]{x1D700}})$ and (b,d,f) logarithm of the temperature fluctuation dissipation rate normalized by its mean value $\log _{10}(\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}/\overline{\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}})$ . For $\unicode[STIX]{x1D6FE}=570$ and $C=0.19$ , at three Stokes numbers: $St=0.003$ (a,b), $St=0.352$ (c,d) and $St=7.343$ (e,f).

We now investigate the effect of the Stokes number, the non-dimensional particle number density, $C$ and the non-dimensional box size, $\unicode[STIX]{x1D6FE}$ , on average quantities characterizing the flow. In figures 3 and 4 we consider the mean turbulent kinetic energy, mean dissipation rate of turbulent kinetic energy, fluid temperature variance and mean rate of dissipation of the temperature fluctuations. It is important to note that in these figures all the quantities are normalized by $\ell _{\ast }$ , $t_{\ast }$ and $\unicode[STIX]{x1D703}_{\ast }$ . Those scales being independent of $St$ , $C$ and $\unicode[STIX]{x1D6FE}$ , it is straightforward to discuss the effect of those parameters.

First we focus in figure 3 on the effect of both Stokes number and mean particle concentration. All quantities appear to peak around $St=1$ , independently of the concentration. For small Stokes numbers, the statistics are dependent on the overall particle number density while, above $St\approx 0.1$ , this dependency vanishes in most cases (the dependency of $\overline{\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}}$ on the concentration is attributed to the small number of particles in the system for $C=0.035$ and $C=0.19$ ).

The influence of the mean particle concentration at low Stokes numbers reflects that the mean inter-particle distance is an important length scale, and that the dynamics is essentially driven by the thermal plumes issued from individual particles. Conversely, for high enough Stokes numbers, the dynamics is not governed by individual particle effects but by particle clustering.

Figure 3.

Influence of Stokes number and the mean particle concentration on the turbulent kinetic energy (a), the dissipation rate of turbulent kinetic energy (b), the fluid temperature variance (c), and the mean dissipation rate of temperature fluctuations (d) normalized respectively by $u_{\ast }=\ell _{\ast }/t_{\ast }$ , $\unicode[STIX]{x1D700}_{\ast }=\ell _{\ast }^{2}t_{\ast }^{-3}$ , $\unicode[STIX]{x1D703}_{\ast }=\unicode[STIX]{x1D6FD}t_{\ast }$ and $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703},\ast }=\unicode[STIX]{x1D703}_{\ast }^{2}t_{\ast }^{-1}$ , in log–log scales. For $\unicode[STIX]{x1D6FE}=83$ and $C=0.035$ (red triangles), 0.19 (black squares), 0.35 (blue circles), 1.82 (green squares), and 8.77 (purple stars).

Figure 4 shows the effect of the Stokes number on the various statistics for five non-dimensionalized box sizes, $\unicode[STIX]{x1D6FE}$ . We observe a maximum at a Stokes number which shifts to higher values when $\unicode[STIX]{x1D6FE}$ is increased, but is in general of order 1. At vanishingly small Stokes numbers, the averaged quantities are practically independent of $\unicode[STIX]{x1D6FE}$ , except for the turbulent kinetic energy. When the Stokes number is increased, the influence of $\unicode[STIX]{x1D6FE}$ is significant and indicates an increase in both velocity and temperature fluctuations with the increasing domain size. This is consistent with the picture of a cluster-dominated system, in which the clusters present wide range of size along with regions depleted in particles of dimensions commensurate with the large scales of the flow. Interestingly though, the domain size has no sizable effect on the mean temperature dissipation rate. Since $\overline{\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}}$ is also roughly independent of the particle concentration, for sufficiently large Stokes number (figure 3), it appears to be only dependent on $St$ when the particles organize in clusters. This has important consequences for the model linking the temperature dissipation rate to the local particle clustering, which we present in the next section.

Figure 4.

Influence of Stokes number and the non-dimensional box size $\unicode[STIX]{x1D6FE}$ on the turbulent kinetic energy (a), the dissipation rate of turbulent kinetic energy (b), the fluid temperature variance (c), and the mean dissipation rate of temperature fluctuations (d) normalized respectively by $u_{\ast }=\ell _{\ast }/t_{\ast }$ , $\unicode[STIX]{x1D700}_{\ast }=\ell _{\ast }^{2}t_{\ast }^{-3}$ , $\unicode[STIX]{x1D703}_{\ast }=\unicode[STIX]{x1D6FD}t_{\ast }$ and $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703},\ast }=\unicode[STIX]{x1D703}_{\ast }^{2}t_{\ast }^{-1}$ , in log–log scales. For $C=0.19$ and $\unicode[STIX]{x1D6FE}=48$ (red triangles), 83 (black squares), 190 (blue circles), 330 (green squares) and 570 (purple stars).

Although the usual definitions of turbulent scales rely on the concept of inertial range, which is not necessarily the case here as discussed in § 4.3, to discuss the turbulent nature of this flow, we provide estimation of the Reynolds number based on the Taylor scale $R_{\unicode[STIX]{x1D706}}=\sqrt{\overline{{u^{\prime }}^{2}}}\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$ . Using the classical estimation of the Taylor scale, $\unicode[STIX]{x1D706}=\sqrt{15\overline{{u^{\prime }}^{2}}\unicode[STIX]{x1D708}/\overline{\unicode[STIX]{x1D700}}}$ , with the characteristic velocity estimated from the turbulent kinetic energy, $\overline{{u^{\prime }}^{2}}\approx 2/3K$ , one obtains the following relation for the Reynolds number based on the Taylor scale: $R_{\unicode[STIX]{x1D706}}=(2/3)\sqrt{15}(K/u_{\ast }^{2})(\overline{\unicode[STIX]{x1D700}}/\unicode[STIX]{x1D700}_{\ast })^{-1/2}$ . Depending on the parameters, $R_{\unicode[STIX]{x1D706}}$ ranges from 15 ( $\unicode[STIX]{x1D6FE}=83$ , $C=8.77$ , $St=0.003$ ) to 250 ( $\unicode[STIX]{x1D6FE}=570$ , $C=0.19$ , $St=1.064$ ) as can be computed from figures 3 and 4. Similarly, the usual Stokes number based on the Kolmogorov time scale, $St_{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ , can be deduced from the value of $\overline{\unicode[STIX]{x1D700}}/\unicode[STIX]{x1D700}_{\ast }$ : $St_{\unicode[STIX]{x1D702}}=St(\overline{\unicode[STIX]{x1D700}}/\unicode[STIX]{x1D700}_{\ast })^{1/2}$ .

4.2 Scaling

Similarly to Hinch (Reference Hinch, Guyon, Nadal and Pomeau1988), Brenner (Reference Brenner1999) and Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002), we consider both thermal and momentum budgets on a volume of fluid containing a cluster of particles in order to represent the active source of heat and buoyancy. Such volumes of fluid are denoted below as ‘effective particles’. The thermal budget of such an effective particle is written as:

where $V_{c}$ is the characteristic volume of the $c^{th}$ cluster, $N_{p,c}$ is the number of actual particles contained in it, $\unicode[STIX]{x1D703}_{c}$ is its characteristic temperature and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D703},c}$ represents its temperature mixing time. The left-hand side of (4.1) is the cluster thermal inertia (recall that here we neglect the thermal inertia of the individual particles), the first term on the right-hand side accounts for the heat source in it and the second term represents the relaxation to thermal equilibrium with the bulk flow. Note that we assume $\overline{\unicode[STIX]{x1D703}}=0$ . i.e. we take the characteristic bulk temperature to be equal to the mean temperature, due to the relatively small particle volume fraction in the cluster. Assuming stationarity, equation (4.1) simplifies to:

The momentum balance for an effective particle reads:

where $u_{c}$ is the velocity of the effective particle and $\unicode[STIX]{x1D70F}_{c}$ is its momentum relaxation time. Notice that the linearized state equation has been used to express the buoyancy term: $g(1-\unicode[STIX]{x1D70C}_{0}/\unicode[STIX]{x1D70C}_{c})=\unicode[STIX]{x1D6FC}g\unicode[STIX]{x1D703}_{c}$ . In case of stationarity one obtains:

Figure 5.

Evolution of statistical quantities with the non-dimensional mean particle concentration $C=\overline{n}\ell _{\ast }^{3}$ for different Stokes numbers, and $\unicode[STIX]{x1D6FE}=83$ . Plots are in log–log scales. (a) the mean particle temperature and comparison with the curve $0.13~C^{-1}$ (dashed lines); (b) the Bolgiano length scale, and the velocity integral length scale in dashed lines; (c) the mean dissipation rate normalized by $\unicode[STIX]{x1D700}_{\ast }=\ell _{\ast }^{2}t_{\ast }^{-3}$ and comparison with the curve $0.43~C^{-3/5}$ (dashed lines); (d) the fluid temperature variance and comparison with the curve $0.2~C^{-4/5}$ (dashed lines); (e) the fluid kinetic energy and comparison with the curve $4~C^{-2/5}$ (dashed lines); (f) mean particle vertical velocity and comparison with the curve $0.07~C^{-4/5}$ (dashed lines).

4.3 Spectral analysis

We present the spectral density of the fluid temperature fluctuations, the turbulent kinetic energy and the particle concentration field in figure 13 for $\unicode[STIX]{x1D6FE}=570$ and $C=0.19$ , which are representative of all domain sizes and concentrations. We consider here the three-dimensional spectra obtained by integration of the Fourier coefficient regardless of the angular dependence. For the temperature fluctuations it reads: $E_{\unicode[STIX]{x1D703}}(k)=\int _{|\boldsymbol{k}|=k}\overline{\unicode[STIX]{x1D703}(\boldsymbol{k})\unicode[STIX]{x1D703}^{\ast }(\boldsymbol{k})}\,\text{d}\boldsymbol{k}$ , with $\unicode[STIX]{x1D703}^{\ast }(\boldsymbol{k})$ the complex conjugate of $\unicode[STIX]{x1D703}(\boldsymbol{k})$ . Also note that here the $\overline{\bullet }$ operator stands for the temporal average. In figure 13 the spectra are normalized by the scale introduced in (2.26)–(2.27) and the wavelength $k$ is normalized by $k_{\ast }=2\unicode[STIX]{x03C0}/\ell _{\ast }$ . A clear power law is present up to wavelengths smaller than $0.1k_{\ast }$ regardless of the Stokes number. For higher wavenumbers, an exponential decay characteristic of the dissipative region is observed. The Stokes number has a dramatic effect on the temperature spectra: at very small particle inertia, the power-law exponent is positive, indicating that the structures of the fields are dominated by the small scales. With increasing Stokes number, the power-law exponent at small wavenumber decreases, going through zero at approximately $St=0.074$ , and reaching an asymptotic value for the largest particle inertia. The decrease of the power-law exponent with $St$ is consistent with the temperature snapshot in figure 1, where we observe that larger structures develop in the temperature field when the Stokes number rises. The evolution of the power-law exponent with the Stokes number is reported in the insert of figure 13. For very small values of the Stokes number, the exponent of the temperature spectra remains constant at approximately $+1/3$ . For $St$ above the threshold, clusters are formed, and then, when the Stokes number is further increased, the clustering of the system is reinforced. This gradual evolution of the particle clustering cause the smooth evolution the spectra exponent, which decreases down to the ‘Obukhov–Corrsin’ value of $-5/3$ obtained for passive scalar fluctuations forced at large scale in isotropic turbulence (Obukhov Reference Obukhov1949; Corrsin Reference Corrsin1951). The continuous evolution of the power-law exponent suggests that it may take non-rational values at intermediate Stokes number. The anomalous values of the exponent are associated with the forcing imposed by the particle clusters, which are seen to have a fractal geometry (as seen in figures 6, 7 and 8). The spectra slope becomes sensitive to the Stokes number from $St=0.019$ at $C=0.19$ which gives an estimate of $St_{0}=0.0036$ consistently with the value proposed for $St_{0,crit}$ in § 4.2.2.

Concerning the turbulent velocity spectra, a somewhat similar behaviour is observed, with a continuous decrease of the power-law exponent from approximately $-1$ at small Stokes numbers to approximately $-5/3$ at large Stokes numbers. The important difference is that the exponent in this case is negative for the entire $St$ range, which implies that the integral $K=\int _{k_{min}}^{k_{max}}E_{u}(k)\,\text{d}k$ does not converge when increasing the size of the domain (see figure 4). On the other hand, the positive power-law exponent, at small Stokes number, for the temperature spectra implies the convergence of the integral $\overline{\unicode[STIX]{x1D703}^{2}}=\int _{k_{min}}^{k_{max}}E_{\unicode[STIX]{x1D703}}(k)\,\text{d}k$ consistent with the independency of the temperature variance on $\unicode[STIX]{x1D6FE}$ (see figure 4) at vanishingly small Stokes number.

The spectral density of the particle concentration field characterizes the scale of the concentration fluctuations which is closely related to the particle clustering. As expected, at small Stokes numbers we observe a $k^{2}$ power law corresponding to a random homogenous particle distribution. For intermediate Stokes number, one observe an almost flat spectrum, indicating that the particle concentration fluctuates at all scales. At large Stokes numbers, most of the concentration fluctuations take place at very large scales while the smallest scales of the concentration fields become randomly distributed (as shown by the $k^{2}$ slope). This evolution of the concentration spectra with $St$ appears to be consistent with the findings of the § 4.2.2.

4.4 Probability density functions

In order to study the fluctuations of the flow, we consider the p.d.f. of some key quantities. The p.d.f. of fluid and particles temperature are given in figure 18. First, we remark that the distributions for both particle temperature and fluid temperature are very close, which is in part a consequence of the assumption of vanishing particle specific heat. At small Stokes number, the temperature p.d.f. deviates somewhat from the Gaussian distribution. Indeed, for vanishing particle inertia, the flow is dominated by the dynamics of individual particles (see § 4.2.1), and as a consequence, the ensemble of discrete hot spots leads to positive fluctuations occurring more frequently than in a normal (Gaussian) process. However, note that with increasing mean particle concentration the temperature fluctuations return to the normal distribution (not shown for brevity). Accordingly, the deviation from normality is interpreted as an effect of the finite number of particles.

For the largest Stokes number, the temperature distribution is closer to Gaussian irrespective of the concentration. This is a consequence of the absence of correlations between the heat sources provided by the ballistic-like motion of the particles with large inertia. For intermediate Stokes numbers, the p.d.f. are strongly skewed, with an exponential tail for the positive values. Similarly, Chillá et al. (Reference Chillá, Ciliberto, Innocenti and Pampaloni1993) observed in a Rayleigh–Bénard flow a transition for the temperature p.d.f. from Gaussian to exponential for Rayleigh number $Ra>10^{7}$ . In our flow this is likely caused by the particle clustering leading to the formation of high temperature regions. This is confirmed by figure 1 where we observe very localized regions of intense positive temperature fluctuations that are well correlated with the high concentration regions.

Figure 18 also presents the p.d.f. of the temperature dissipation rate $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}$ . For small Stokes numbers, the fluctuations of $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}$ remain rather small reflecting the fact that the temperature gradient fields are relatively ‘smooth’, as seen in figure 2. When the Stokes number is increased, extreme fluctuations of $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}$ develop and an increasingly stretched tail of the p.d.f. is observed. This is consistent with the very intense scalar dissipation rate organized on thin and long line, also visible in figure 2. Those lines have a thickness of the order of the dissipative length scale and length of the order of the integral scale. The very large values of $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703}}$ observed at large Stokes numbers imply that a significant part of the scalar dissipation occurs in those very concentrated regions, which is the manifestation of the small-scale intermittency. Finally, for the largest Stokes number the fluctuations of the temperature dissipation rate are reduced, as the strongly inertial particles are less preferentially concentrated and temperature field is smoothed out by the large-scale mixing.

Before focusing on the velocity p.d.f., we present in figure 19 the evolution of the variance of both horizontal and vertical velocity components with the Stokes number for the various mean particle concentration. In this plot both velocity components are normalized by the turbulent kinetic energy $K$ . It is observed that the vertical velocity has a greater variance than the horizontal component. Basically we note that the ratio between the components is approximately 3, although the difference is slightly more pronounced at small Stokes numbers. Furthermore, we observe only a small dependence of this ratio on the average particle concentration. Moreover, we note that the effect of $\unicode[STIX]{x1D6FE}$ is very small (not shown).

Figure 19.

Evolution with the Stokes number of the variance of the horizontal (continuous lines) and vertical (dashed lines) velocity components both normalized by the mean turbulent kinetic energy $K$ . For $\unicode[STIX]{x1D6FE}=83$ and $C=0.035$ (red triangles), $0.19$ (black squares), $0.35$ (blue circles), $1.82$ (green squares), $8.77$ (purple stars).

The p.d.f. of the fluid velocity is given in figure 20. We consider both horizontal and vertical velocity components normalized by their standard deviation. We observe that, although the highest anisotropy of the turbulent structures is seen at low $St$ (figures 17 and 19), the p.d.f. of horizontal and vertical velocities are actually both very close to the normal distribution. This suggests that at low $St$ the anisotropy of the velocity vector components is mainly due to their difference in variance. On the other hand, for larger $St$ , only the vertical component remains close to the normal distribution whereas the horizontal component present much larger fluctuations. This indicates that the large-scale anisotropy of the structure depicted in figure 17, produces large-scale intermittency in addition to the difference between the horizontal and vertical velocity variances.

Figure 20.

(a) P.d.f. of the fluid velocity: horizontal component (continuous lines) and vertical component (dotted line). For the different Stokes numbers ( $St=0.003$ , 0.019, 0.074, 0.352, 1.064, 7.343 and 29.36) each shifted upward by two decades respectively. (b) P.d.f. of the dissipation (continuous lines) and of the production (dotted line), for the different Stokes numbers each shifted upward by two decades respectively. With $\unicode[STIX]{x1D6FE}=570$ and $C=0.19$ . Comparison with the Gaussian distribution in grey dashed lines.

Figure 20 also reports the comparison between the p.d.f. of the dissipation, $\unicode[STIX]{x1D700}$ , and of the production of the turbulent kinetic energy, $\unicode[STIX]{x1D6FC}g\unicode[STIX]{x1D703}w$ , both normalized by the mean dissipation. We observe that when the Stokes number is increased the dissipation p.d.f. develop very stretched tails. Similarly to the temperature dissipation rate, the very large fluctuations of $\unicode[STIX]{x1D700}$ are connected with the intermittent nature of the dissipation process at large Stokes number, as seen in the snapshots of the figure 2. Regarding the production term, the occurrence of very large fluctuations increases up to a Stokes number of order $St=O(1)$ , and then the tails become narrower when the Stokes number is further increased. This reflects the nature of the forcing, due to both the discrete particle effect and the very strong correlation between the high temperature region and the large raising velocity. When the Stokes number is high enough to have very important clustering, this correlation is reinforced, and the forcing is mainly operated by the few largest clusters. The particles of greater inertia can instead escape from the regions of high temperature leading to a more spatially homogeneous forcing. Also we remark that the production term, although positive on average, also takes large negative values, corresponding to reverse transfer of turbulent kinetic energy into potential energy.

5.1 Effect of finite particle thermal capacity and momentum two-way coupling

In the previous sections the transfer of momentum from the particle phase to the fluid phase has been neglected, an assumption which is typically made for systems of very small particles and for limited number density. Likewise, the particles have been so far considered in thermal equilibrium with the surrounding fluid by assuming zero thermal inertia ( $c_{p}=0$ ). This is also a reasonable assumption when extremely small particles are considered. In this section we investigate the effects of the finite thermal inertia of the particles along with the momentum feedback from the particles to the fluid. For constant particle number density, both effects are promoted when the Stokes number is increased, as seen from (2.28).

In order to discuss the influence of the two aforementioned effects, we consider therefore two additional cases. In the first one we fix the solid-to-fluid specific heat ratio as $c_{p}/c_{f}=4.18$ but we disregard the effect of the particle momentum feedback on the fluid phase. In the second case we keep the finite particle thermal inertia to $c_{p}/c_{f}=4.18$ and in addition we set full momentum and heat transfer coupling. For both cases the parameters are $C=0.35$ , $\unicode[STIX]{x1D6FE}=80$ , $\unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{0}=909$ and $Pr=1$ , and the same range of Stokes numbers as in the previous sections is considered (from $St=0.003$ to $St=29.36$ ) as summarized in table 2. These parameters correspond to a mass fraction ranging from $\unicode[STIX]{x1D6FC}_{m}=7.5\times 10^{-5}$ to $\unicode[STIX]{x1D6FC}_{m}=9.88\times 10^{-1}$ . These two sets of simulations are then compared to the simplified case considered in the previous section. In all these cases we impose the same value of $\unicode[STIX]{x1D6FD}$ , nevertheless it should be noted that when $c_{p}$ is finite, the energy injected in the system increase substantially with the mass fraction as seen from (2.12) since most of the energy input is required to heat the dispersed phase. The total energy input scales as $\unicode[STIX]{x1D6FD}(1+\unicode[STIX]{x1D712})$ . As apparent from the definition of $St$ , for $St=29.36$ the particle diameter is $d_{p}/\ell _{\ast }=\sqrt{18St\unicode[STIX]{x1D70C}_{0}/\unicode[STIX]{x1D70C}_{p}}\approx 0.76$ which make the assumption $d_{p}\ll \ell _{\ast }$ somehow questionable. However, even the most inertial particles would still be fairly small compared to the energy containing structures, and therefore size effects are not expected to qualitatively change our findings. Also note that for the more realistic case (i.e. finite particle thermal inertia and momentum two-way coupling) we considered two more intermediate values for the Stokes number: $St=4.26$ and $St=12.69$ .

Figure 21 shows the evolution of the turbulent kinetic energy, the variance of the fluid temperature and the energy and temperature dissipation rate with the Stokes number, for the three cases discussed above (with and without particle thermal inertia and with and without momentum two-way coupling). As expected, at small Stokes numbers, these statistical quantities are identical in the three cases since the mass fraction is vanishingly small. At larger Stokes number, when only the finite thermal inertia of the particles is accounted for, the temperature field appears more homogeneous as is visible from the decrease of both temperature variance and mean temperature dissipation rate presented in figure 21. Indeed, due to the finite thermal inertia of the particles, the instantaneous total inter-phase heat exchange presents temporal fluctuations, whereas neglecting the thermal inertia of the particles results in a constant heat release in the fluid, as seen in (2.14). And since for large particle thermal inertia the particle temperature could be considered as quasi-constant, the inter-phase heat transfer is promoted when the particles are in regions where the fluid is cold, and is decreased in hot regions. This effect tends to smooth the temperature fields. When the Stokes number is of order 1, as previously observed, the particles present a very intense clustering and the residence time of particles in the large clusters is much longer than their thermal relaxation time. As a consequence, the particles can relax towards thermal equilibrium and the mechanical behaviour of the system is only slightly affected by the thermal inertia of the particles, as is visible from the plot of the turbulent kinetic energy and turbulent dissipation rate.

Figure 21.

Evolution with the Stokes number of the turbulent kinetic energy (a), the dissipation rate of turbulent kinetic energy (b), the fluid temperature variance (c), and the mean dissipation rate of temperature fluctuations (d) normalized respectively by $u_{\ast }=\ell _{\ast }/t_{\ast }$ , $\unicode[STIX]{x1D700}_{\ast }=\ell _{\ast }^{2}t_{\ast }^{-3}$ , $\unicode[STIX]{x1D703}_{\ast }=\unicode[STIX]{x1D6FD}t_{\ast }$ and $\unicode[STIX]{x1D700}_{\unicode[STIX]{x1D703},\ast }=\unicode[STIX]{x1D703}_{\ast }^{2}t_{\ast }^{-1}$ . At $\unicode[STIX]{x1D6FE}=83$ and $C=0.35$ and for (i) $c_{p}/c_{f}=0$ without momentum two-way coupling (red triangles); (ii) $c_{p}/c_{f}=4.18$ without momentum two-way coupling (blue circles) and (iii) $c_{p}/c_{f}=4.18$ with momentum two-way coupling (black squares).

In contrast, the attenuation of the turbulent kinetic energy and turbulent dissipation rate is more pronounced at $St=O(1)$ when accounting for the momentum two-way coupling. This is explained by the fact that the large inertia of the clusters tends to drag the fluid. For larger Stokes number ( $St>1$ ), that is to say for large mass fractions, the case with momentum two-way coupling presents a sudden increase of both $K$ and $\overline{\unicode[STIX]{x1D700}}$ . This second peak is attributed to the forcing in the fluid momentum equation operated by the particle feedback. At that Stokes number the mass fraction becomes important and the ballistic motion of the particles causes the fluid to be dragged by the particles, providing an efficient source of turbulent agitations as well as temperature mixing. We would like to draw attention to the fact that in order to assess the influence of relaxing the assumptions explored here, quantities in figure 21 are normalized by $\ell _{\ast }$ , $t_{\ast }$ and $\unicode[STIX]{x1D703}_{\ast }$ , which imply implicitly that $\unicode[STIX]{x1D6FD}$ is kept fix between all of the explored cases. However, as stated above, a more meaningful comparison could be to consider a constant total energy input, i.e. $\unicode[STIX]{x1D6FD}(1+\unicode[STIX]{x1D712})=\text{cst}$ .

Figure 22.

Snapshots in a $x-z$ plane: (a) of the logarithm of the particle concentration field (as computed by (3.1)) normalized by the mean concentration $n/\overline{n}$ , and (b) the fluid temperature fluctuations normalized by its standard deviation $\unicode[STIX]{x1D703}/\overline{\unicode[STIX]{x1D703}^{2}}^{1/2}$ . For $\unicode[STIX]{x1D6FE}=80$ and $C=0.35$ , $St=0.352$ and $Fr=0.09$ .

Finally, we remark that when the mass fraction is small enough to neglect the inter-phase momentum two-way coupling, one can also neglect the finite thermal inertia of the particle. Of course this conclusion only holds for not too high solid-to-fluid specific heat ratio.

5.2 Effect of the finite particle Froude number

Gravity has so far been considered previously only as a source of buoyancy, without directly affecting the particles, i.e. $Fr=\infty$ in (2.34). We here consider cases for various Froude numbers ( $Fr=0.09$ , $0.32$ , $0.48$ , $3.16$ , $10$ and $\infty$ ), keeping all other parameters constant. In particular, we impose $\unicode[STIX]{x1D6FE}=80$ , $St=0.3$ , $C=0.35$ $Pr=1$ and $\displaystyle \unicode[STIX]{x1D70C}_{p}/\unicode[STIX]{x1D70C}_{0}=909$ , $\displaystyle c_{p}/c_{f}=0$ . The parameters of this set of simulations are given in table 2.

Snapshots of the particle local concentration and fluid temperature corresponding to $Fr=0.09$ are presented in figure 22. It is observed that for this small Froude number, the particles are organized in very long vertical streaks. It is to note that similar patterns have been observed experimentally in sedimentation of both positively and negatively buoyant particles (Weiland, Fessas & Ramarao Reference Weiland, Fessas and Ramarao1984; Davis & Acrivos Reference Davis and Acrivos1985). This peculiar particle distribution significantly departs from the $1/Fr=0$ case in figures 1. The elongated particle streaks are associated with elongated plumes of high temperature fluid.

Figure 23 shows the evolution of the mean vertical particle velocity with the Froude number. The particle velocity decreases linearly with $1/Fr^{2}$ , and it is negative for $Fr<0.7$ . On the other hand, for $Fr>1$ the (positive) mean vertical particle velocity is independent of the Froude number, confirming that the limit $Fr\rightarrow \infty$ considered previously is meaningful. The fact that the transition between falling and rising particles occurs at a Froude number of order 1 further corroborates the relevance of the characteristic scales introduced in (2.25)–(2.27). We remark that, in the regime in which particles settle ( $Fr<0.7$ ), their fall speed is somewhat smaller than the Stokesian settling velocity, $u_{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D70F}_{p}g$ . This is explained by the fact that particles cluster and create hot plumes of upward moving fluid, a phenomenon also observed in settling of hot particles through forced turbulence (Frankel et al. Reference Frankel, Pouransari, Coletti and Mani2016).

In figure 23 the particle temperature versus the Froude number is also plotted. Note that, since $c_{p}=0$ for this set of simulations, this corresponds to the fluid temperature at the particle positions ( $\unicode[STIX]{x1D703}_{p}=\unicode[STIX]{x1D703}(\boldsymbol{x}_{p})$ ). Similarly to the mean settling velocity, the particle temperature is not strongly influenced by the Froude number when $Fr>1$ . For small Froude numbers the particle temperature decreases dramatically, due to the enhanced mixing introduced by the particle settling, but on average it remains higher than the fluid temperature ( $\langle \unicode[STIX]{x1D703}\rangle _{p}>0$ ). The local excess of temperature at the particle location leads to the upward, buoyancy-driven motions mentioned above, which hinder the particle settling. The small peak observed in the particle temperature around $Fr\approx 1$ can be attributed to the fact that at this regime the mean vertical velocity of the particles is almost zero, leading to higher residence time of the particles in the plumes. In a recent study of DNS of hydrodynamically forced turbulence laden with heated inertial particles (Frankel et al. Reference Frankel, Pouransari, Coletti and Mani2016), it was shown that the buoyant plume shed by the clusters breaks down the preferential sweeping mechanism by which particle settle faster as they favour the downward side of turbulent eddies (Wang & Maxey Reference Wang and Maxey1993).

Figure 23.

Evolution with the Froude number of (a) the mean vertical particle velocity and comparison with the Stokes terminal velocity $u_{\unicode[STIX]{x1D70F}}=\unicode[STIX]{x1D70F}_{p}g$ (dashed lines) and (b) the mean particle temperature, for $\unicode[STIX]{x1D6FE}=80$ and $St=0.3$ $C=0.35$ , normalized by the critical values $u_{\ast }^{2}$ and $\unicode[STIX]{x1D703}_{\ast }$ . Inserts: evolution with $Fr^{-2}$ in linear scales.

For the small Froude number cases in which very elongated structures appear, the present domain size in vertical direction may not be sufficient to prevent biases due to the periodic boundary conditions. Further studies of these relevant regimes are warranted.