1 Introduction
Multiphase flows in which a dense dispersed phase interacts with a lighter carrier fluid are ubiquitous in nature and industrial applications. Often the fluid flow is in the turbulent regime, and the inertial particles cannot follow its rapid velocity fluctuations. This velocity lagging can lead to high local concentration in zones of high strain and away from vortex cores (Maxey Reference Maxey1987; Lazaro & Lasheras Reference Lazaro and Lasheras1989; Squires & Eaton Reference Squires and Eaton1991). This phenomenon, known as preferential concentration or turbulence clustering, has been widely investigated and is considered a critical mechanism for the growth of cloud droplets (Shaw Reference Shaw2003; Grabowski & Wang Reference Grabowski and Wang2013), the dust agglomeration in circumstellar nebulas (Cuzzi et al. Reference Cuzzi, Hogan, Paque and Dobrovolskis2001) and the mixing of air and fuel in spray combustion (Sahu, Hardalupas & Taylor Reference Sahu, Hardalupas and Taylor2014).
For particles much denser than the carrier fluid and small compared to the flow scales, the relevant parameter describing the particle–fluid interaction is the Stokes number
$St={\it\tau}_{p}/{\it\tau}_{f}$
, the ratio between the Stokes aerodynamic response time of the particle
${\it\tau}_{p}={\it\rho}_{p}D_{p}^{2}/18{\it\mu}$
and the fluid time scale
${\it\tau}_{f}$
. When the latter is taken as the Kolmogorov time scale
${\it\tau}_{{\it\eta}}$
, the maximum level of preferential concentration is typically found at
$St_{{\it\eta}}\approx 1$
(Wang & Maxey Reference Wang and Maxey1993; Eaton & Fessler Reference Eaton and Fessler1994), although there is evidence that particles of
$St_{{\it\eta}}\approx 10$
are clustered at larger scales (Yoshimoto & Goto Reference Yoshimoto and Goto2007). Particles of
$St_{{\it\eta}}=1$
also display maximum dispersion rates and maximum time rate of change of their temperature fluctuations (Wetchagarun & Riley Reference Wetchagarun and Riley2010). In the last decade, experimental, numerical and theoretical studies have deepened our understanding of the origins of turbulence clustering (Goto & Vassilicos Reference Goto and Vassilicos2008; Salazar et al.
Reference Salazar, De Jong, Cao, Woodward, Meng and Collins2008; Gibert, Xu & Bodenschatz Reference Gibert, Xu and Bodenschatz2012; Bragg & Collins Reference Bragg and Collins2014), its consequences on the statistics of particle motion (Ayyalasomayajula et al.
Reference Ayyalasomayajula, Gylfason, Collins, Bodenschatz and Warhaft2006; Bec et al.
Reference Bec, Biferale, Boffetta, Celani, Cencini, Lanotte, Musacchio and Toschi2006) and on the particle collision rate (Sundaram & Collins Reference Sundaram and Collins1997; Wang, Wexler & Zhou Reference Wang, Wexler and Zhou2000).
The dynamics of gravitational settling of inertial particles through turbulence has received less attention compared to turbulence clustering despite its strong practical relevance. Turbulence can increase the settling velocity of small inertial particles beyond the nominal terminal velocity
$V_{t}={\it\tau}_{p}g$
in a quiescent fluid because trajectories that sample the downward side of the eddies are favoured. This mechanism, referred to as preferential sweeping, was first demonstrated in homogeneous isotropic turbulence by Wang & Maxey (Reference Wang and Maxey1993). Measurements from Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002) and Yang & Shy (Reference Yang and Shy2003, Reference Yang and Shy2005) confirmed this effect, which is especially strong when
$St_{{\it\eta}}\approx 1$
and
$V_{t}/u_{rms}\approx 0.5$
,
$u_{rms}$
being the root-mean-square (r.m.s.) fluid velocity fluctuation (Yang & Lei Reference Yang and Lei1998; Dejoan & Monchaux Reference Dejoan and Monchaux2013). Depending on which turbulent scales are considered crucial for the settling rate, researchers have used the ratios
$V_{t}/u_{rms}$
or
$V_{t}/u_{{\it\eta}}$
(Yang & Lei Reference Yang and Lei1998; Yang & Shy Reference Yang and Shy2003). For an in depth discussion on the non-dimensional parameters governing the settling of small particles in vortical flows, the reader is referred to Davila & Hunt (Reference Davila and Hunt2001).
Although they are less commonly described in the literature, mechanisms exist by which turbulence can instead reduce the particle settling velocity. Very small, tracer-like particles can be trapped in vortices (Tooby, Wick & Isaacs Reference Tooby, Wick and Isaacs1977), while large particles can be slowed down by the nonlinear drag due to fluid velocity fluctuations (Mei Reference Mei1994). Moreover, heavy particles falling rapidly through turbulence might spend relatively more time in the upward-moving sides of the vortices, which also hinders their settling, a phenomenon referred to as loitering (Nielsen Reference Nielsen1993; Good, Gerashchenko & Warhaft Reference Good, Gerashchenko and Warhaft2012). Recent results of Good et al. (Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014) indicate that settling velocity reductions for heavy particles in direct numerical simulations are only possible when a nonlinear drag correction is introduced.
The above studies considered only the effects of the flow on the particles, neglecting the reaction force on the fluid. If the loading is sufficiently high or the particle size is not small compared to the smallest flow scales, a two-way coupled regime is established. Most studies concerned with two-way coupling consider zero-gravity conditions, and find that turbulence is either enhanced or reduced depending on the particle size, Stokes number, loading and settling rate (Yang & Shy Reference Yang and Shy2005; Poelma, Westerweel & Ooms Reference Poelma, Westerweel and Ooms2007; Balachandar & Eaton Reference Balachandar and Eaton2010). For particles of diameter
$D_{p}$
much smaller than the Kolmogorov scale
${\it\eta}$
and Stokes number
$St_{{\it\eta}}\approx 1$
, studies in both homogenous and sheared turbulent flows are consistent in showing a global decrease of turbulence, with depletion of turbulent energy at large scales and injection at small scales, an effect referred to as ‘pivoting’ or ‘cross-over’ (Elghobashi & Truesdell Reference Elghobashi and Truesdell1993; Boivin, Simonin & Squires Reference Boivin, Simonin and Squires1998; Poelma et al.
Reference Poelma, Westerweel and Ooms2007; Gao, Li & Wang Reference Gao, Li and Wang2013; Gualtieri et al.
Reference Gualtieri, Picano, Sardina and Casciola2013).
When gravity acts on a two-way coupled dispersed flow, the breaking of symmetry causes anisotropy in otherwise isotropic environments. Heavy particles settling through turbulence may exert a non-negligible collective drag on the carrier fluid, which adds to the preferential sweeping and increases the settling rate (Bosse, Kleiser & Meiburg Reference Bosse, L. and Meiburg2006), while the turbulent energy is redistributed from the horizontal to the vertical direction (Ferrante & Elghobashi Reference Ferrante and Elghobashi2003; Poelma et al. Reference Poelma, Westerweel and Ooms2007). In the absence of hydrodynamic turbulence forcing, fluid velocity fluctuations will be caused solely by the momentum coupling with the settling particles. The induced eddies can span a wide spectrum of scales due to the formation of large clusters (Kajishima & Takiguchi Reference Kajishima and Takiguchi2002; Capecelatro, Desjardins & Fox Reference Capecelatro, Desjardins and Fox2014). Similarly, rising bubbles induce anisotropic turbulence in the liquid phase (so-called pseudo-turbulence), with higher energy content in the vertical direction (Lance & Bataille Reference Lance and Bataille1991; Bunner & Tryggvason Reference Bunner and Tryggvason2002; Mazzitelli & Lohse Reference Mazzitelli and Lohse2009). Recently we have shown that, in an otherwise stagnant particle–gas mixture, the radiative heating of particles can result in the formation of buoyant plumes and highly intermittent concentration, leading to a turbulent state sustained only by the radiation (Zamansky et al. Reference Zamansky, Coletti, Massot and Mani2014). We have also shown that the preferential concentration of particles by turbulence can impact the heat transfer between the two phases in such systems (Pouransari, Mortazavi & Mani Reference Pouransari, Mortazavi and Mani2015).
In the present study we explore the behaviour of heated particles settling through homogeneous turbulence and the interaction between dispersed and continuous phases. The physical parameters are chosen to be representative of a dilute, optically thin mixture of air and solid particles subject to thermal radiation. The particle Stokes number and settling velocity are chosen to produce significant preferential concentration. Numerical simulations are performed taking into account the two-way coupling of both momentum and temperature between the dispersed and continuous phase. Since particles much smaller than the smallest flow scales are considered, all relevant spatial and temporal scales can be resolved. A simple schematic is shown in figure 1. To the best of our knowledge, this regime has not been investigated before.
We consider the scenario in which small particles are settling under the influence of gravity and are externally irradiated (a). When they are heated, the clustered particles shed buoyant plumes which slow down their descent (b).

2 Methodology
2.1 Governing equations
The variable-density low-Mach Navier–Stokes equations for a fluid of constant viscosity
${\it\mu}$
are
where
$u_{i}$
are the velocity components,
$t$
is time,
$x_{i}$
are the spatial coordinates,
$p$
is the hydrodynamic pressure,
${\it\tau}_{ij}={\it\mu}((\partial u_{i}/\partial x_{j})+(\partial u_{j}/\partial x_{i})-(2/3)(\partial u_{k}/\partial x_{k}){\it\delta}_{ij})$
is the Newtonian viscous stress tensor,
${\it\rho}$
is the fluid density,
$g_{i}$
(
$i=3$
) is the gravitational acceleration,
${\it\delta}$
is the Dirac delta function,
$N_{p}$
is the number of particles and
$v_{i,n}$
and
$y_{i,n}$
are the velocity and spatial coordinates of the
$n$
th particle. Though the domain is triply periodic, the addition of gravity renders the vertical direction asymmetric with respect to the two horizontal directions. Because we consider spherical Stokesian particles, their mass and aerodynamic response time are
$m_{p}={\it\rho}_{p}{\rm\pi}D_{p}^{3}/6$
and
${\it\tau}_{p}={\it\rho}_{p}D_{p}^{2}/18{\it\mu}$
, respectively,
${\it\rho}_{p}$
being the particle density and
$D_{p}$
the particle diameter. The fluid is considered to be an ideal gas of constant thermal conductivity
$k$
and heat capacities
$C_{v}$
and
$C_{p}$
(at constant volume and pressure, respectively), for which the energy equation reads
where
$T_{f}$
is the fluid temperature,
$T_{p,n}$
is the temperature of the
$n$
th particle and
$h$
is the convective heat transfer coefficient, which for a Stokesian particle can be calculated from the Nusselt number
$Nu=hD_{p}/k=2$
. The heat transfer from each particle is evaluated using the fluid temperature in the vicinity of the particle in a manner analogous to the evaluation of the momentum two-way coupling term. The gas properties reflect those of air at ambient temperature and pressure. Particles are individually tracked along their trajectories, according to the simplified particle equation of motion, where only contributions from Stokes drag and gravity are retained,
Under the Stokesian particle assumption, the nominal terminal velocity of the particles is
$V_{t}=g{\it\tau}_{p}$
. Each particle is subject to a radiative heat flux
$I_{o}$
. The carrier phase is transparent to radiation, whereas the incident radiative flux on each particle is completely absorbed. Because we focus on relatively small volume fractions, the fluid–particle medium is considered optically thin. Under these hypotheses, the direction of the radiation is inconsequential, each particle receives the same radiative heat flux, and its temperature
$T_{p}$
is governed by
where
$C_{v,p}$
is the particle specific heat, which is also assumed to be constant with respect to temperature.
Several non-dimensional parameters are introduced in order to characterize the system dynamics. The density ratio
${\it\rho}_{p}/{\it\rho}_{o}$
measures the relative density of the particle bulk material to the density of the fluid, and the mass loading ratio
${\it\phi}_{m}$
measures the total mass of the particles compared to the mass of the fluid. These parameters were adjusted in the simulations by increasing the density of the particle material. The temperature ratio
$T_{p}/T_{o}$
measures the expected temperature of the particle phase compared to the fluid, and was adjusted by increasing the radiation intensity
$I_{o}$
. The Stokes number
$St_{{\it\eta}}$
measures the response time of the particle compared to the Kolmogorov time scale. Increasing the particle density
${\it\rho}_{p}$
increases the particle relaxation time, but the added radiative input also decreases the Kolmogorov time scale, further increasing the value of the Stokes number. Additional details, including the non-dimensional formulation of the governing equations, may be found in appendix A.
2.2 Numerical implementation
The application that originally motivated our study is solid particle solar receivers, a technology for collecting, storing and transforming solar energy (typically through focusing arrays) onto a stream of gas laden with absorbing particles (see for example the review of Tan & Chen (Reference Tan and Chen2010)). As such we investigate the dynamics of the gas and turbulence in regimes typical of these applications, with material densities
${\it\rho}_{p}$
between 1000 and
$10\,000~\text{kg}~\text{m}^{-3}$
and radiation intensity levels
$I_{o}$
between 0 and
$2~\text{MW}~\text{m}^{-2}$
. Although the particles in such applications are irregular, we approximate the particles as spherical and mono-disperse in order to simplify the analysis. We also consider the mixture of gas and particles in a periodic geometry in order to remove the effects of the physical configuration and isolate the relevant physics. Additionally, the regimes considered here may find relevance in high temperature sooting flames, in which small particles are produced and exposed to radiation intensities comparable to the ones studied here.
Equations (2.1)–(2.3) above are solved using a staggered grid formulation and second-order central finite differences. The choice of a finite difference scheme was made in order to easily handle the effects of variable density and particle tracking using projections and interpolations, which would otherwise be complex to handle in a spectral method. Equations (2.4) and (2.5) are solved for each individual particle in the domain. Both the fluid and particle equations are integrated in time using a fourth-order Runge–Kutta time integrator. To compute heat and momentum exchange between the fluid and particle phase, Eulerian temperature and velocity fields are evaluated at the location of the particles using tri-linear interpolation as an approximation to the Dirac delta function in (2.2) and (2.3). Coupling from particles to the fluid is performed with tri-linear projection back to the nodes of the cell containing each particle. The initial conditions for the simulations are generated from a Passot–Pouquet spectrum with the particles distributed randomly throughout the computational domain. The turbulence is maintained using a linear forcing method (Lundgren Reference Lundgren2003; Rosales & Meneveau Reference Rosales and Meneveau2005). The domain is a triply periodic cube of size
$L_{o}=0.237~\text{m}$
, which is uniformly discretized in
$N^{3}$
grid points such that the Kolmogorov scale is resolved in each simulation case (
$L_{o}/N={\it\Delta}\sim {\it\eta}$
).
To enable separate assessment of radiation effects and momentum two-way coupling, we performed additional simulations in which no radiation was added and a simulation in which only one-way momentum coupling between the fluid and particles was allowed. For these simulations we considered a Taylor microscale Reynolds number of
$Re_{{\it\lambda}}=65$
, which can be resolved with
$N=256$
using cubic meshes. For the cases in which two-way coupling enhances the turbulence agitation (and so reduces the Kolmogorov length),
$N$
is increased accordingly. The fluid is considered an ideal gas, with constant dynamic viscosity
${\it\mu}=1.9\times 10^{-5}~\text{Pa}~\text{s}$
and an initial density of
${\it\rho}_{o}=1.2~\text{kg}~\text{m}^{-3}$
at a temperature of
$T_{o}=300~\text{K}$
. The inhomogeneous heat transfer between the dispersed and continuous phases leads to a spatial variation in the gas density. Therefore, the gas aerodynamic pressure is governed by a variable coefficient Poisson equation which is solved via an iterative spectral method (Pares-Sierra & Vallis Reference Pares-Sierra and Vallis1989) at every substep. The mean fluid momentum is set to zero in each time step, which is equivalent to applying a mean hydrostatic pressure gradient to the fluid. This is necessary to prevent the fluid from accelerating continuously in the vertical direction due to the drag imposed by the settling particles and the buoyancy (Maxey & Patel Reference Maxey and Patel2001; Bosse et al.
Reference Bosse, Kleiser, Hartel and Meiburg2005). Similarly, the mean temperature of the fluid is artificially set to be constant to prevent the fluid from warming up indefinitely.
Physical parameters and main flow statistics. For all cases
$D_{p}=40~{\rm\mu}\text{m}$
,
$N_{p}=3.99\times 10^{6}$
,
${\it\Phi}_{v}=10^{-5}$
,
${\it\mu}=1.9\times 10^{-5}~\text{Pa}~\text{s}$
,
$k=0.0275~\text{W}~\text{m}^{-1}~\text{K}^{-1}$
,
$C_{v}=723~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$
,
$C_{p}=1012~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$
, and
$C_{vp}=700~\text{J}~\text{kg}^{-1}~\text{K}^{-1}$
.

2.3 Physical parameters
The main parameters describing the turbulent fluid phase and the dispersed particle phase for the performed simulations are listed in table 1. The values are representative of the system after it has become stationary (typically after ten large-eddy turnover times). Because of the anisotropy caused by the two-way coupling, we define as the velocity and length scale of the large eddies, respectively,
where
$\text{TKE}=(2u_{rms}^{2}+w_{rms}^{2})/2$
is the turbulent kinetic energy,
$u_{rms}$
and
$w_{rms}$
are the r.m.s. fluctuations of the horizontal and vertical components and
${\it\epsilon}$
is the turbulent dissipation rate. The hydrodynamic artificial forcing is the same for all cases, and produces a microscale Reynolds number
$Re_{{\it\lambda},{\it\Phi}_{v}=0}=65$
in the one-way coupled regime (
${\it\Phi}_{v}=0$
subscript denotes no volume loading). In every case the particle diameter is assumed to be
$D_{p}=40~{\rm\mu}\text{m}$
, which results in
$D_{p}/{\it\eta}\lesssim 0.1$
and
$Re_{p}=|v_{p}-u|D/{\it\nu}\lesssim 1$
, confirming the validity of the point-particle approximation. The volume fraction is
${\it\Phi}_{v}=10^{-5}$
(corresponding to a mass fraction
${\it\Phi}_{m}$
between 1 % and 10 %), which is considered in the two-way coupled regime (Elghobashi Reference Elghobashi1994). The momentum transfer from the heavy particles and the buoyancy-driven plumes strongly modify the turbulence structure, as evident for example from the changes in
$Re_{{\it\lambda}}$
between the various cases. The change in turbulent scales also alters the effective Stokes number. The level of radiation is associated with a steady state particle temperature increase with respect to the surrounding fluid, which from (2.5) is
where
$Nu=hD_{p}/k=2$
has been substituted, which is valid for spherical particles at the present
$Re_{p}$
range. The nominal particle-to-fluid temperature ratio is 1.3 for the intermediate and 2.2 for the higher radiation level. The use of the temperature ratio to represent the thermal input to the system is a direct measure of the temperatures present in the flow and also generalizes the configuration to cases where the particles are heated due to other effects besides radiative transfer.
3 Results
3.1 Settling of non-heated particles
3.1.1 Particle clustering
For all considered cases, the particle Stokes number indicates that strong turbulence clustering is expected. This is confirmed by analysing the probability density function (PDF) of the volumes of the three-dimensional Voronoi cells around each particle, which are defined based on the distance to the neighbouring particles. The Voronoi cell size distribution is commonly used to recognize the presence of clusters/voids, i.e. regions of significantly higher/lower concentration with respect to the spatial average over the entire domain (Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2010). Figure 2 shows the PDF of Voronoi cells normalized by the average volume
$\langle V\rangle$
. The PDFs for the inertial particles are much broader than the
${\it\Gamma}$
distribution expected for randomly distributed particles (Ferenc & Neda Reference Ferenc and Neda2007; Tagawa et al.
Reference Tagawa, Mercado, Prakash, Calzavarini, Sun and Lohse2012), which is a hallmark of turbulence clustering. The variance of the Voronoi cell volumes
${\it\sigma}_{V}^{2}$
, which is often used to quantify the degree of clustering (Monchaux et al.
Reference Monchaux, Bourgoin and Cartellier2010; Dejoan & Monchaux Reference Dejoan and Monchaux2013), is also larger for the heavier particles, which is expected given the higher preferential concentration at the intermediate Stokes numbers.
Voronoi volume distributions for the non-radiated particles with different density ratios, full circles:
${\it\rho}_{p}/{\it\rho}_{o}=833$
(
$St_{{\it\eta}}=0.17$
,
${\it\sigma}_{V}=0.86$
), open circles:
${\it\rho}_{p}/{\it\rho}_{o}=4166$
(
$St_{{\it\eta}}=1.6$
,
${\it\sigma}_{V}=2.3$
), triangles:
${\it\rho}_{p}/{\it\rho}_{o}=8333$
(
$St_{{\it\eta}}=5.1$
,
${\it\sigma}_{V}=2.5$
), line:
${\it\Gamma}$
distribution.

3.1.2 Momentum coupling
In cases where the settling speed is much smaller than the fluid velocity fluctuations, the momentum two-way coupling is known to reduce the turbulent kinetic energy. Bosse et al. (Reference Bosse, L. and Meiburg2006) performed two-way coupled simulations in the presence of gravity, and found that the presence of particles enhanced the dissipation and decreased the turbulent kinetic energy. They considered cases in which
$V_{t}/u_{{\it\eta}}$
and
$V_{t}/u_{rms}$
were of order 1 and 0.3, respectively. In contrast, we consider particles of significant settling velocities compared to both the small- and large-scale turbulent velocity fluctuations, and we find that the turbulence intensity is enhanced by the faster settling particles. This is clearly illustrated in figure 3, that displays normalized values of TKE and
${\it\epsilon}$
. For
${\it\rho}_{p}/{\it\rho}_{o}=833$
, which corresponds to settling velocities of the same order as the turbulent velocity scales (see table 1), the turbulent energy and dissipation rate are slightly smaller than those in the one-way coupled flow, in agreement with previous findings (Bosse et al.
Reference Bosse, L. and Meiburg2006). However, for heavier particles with settling parameters
$V_{t}/u_{{\it\eta}}=O(10)$
, both TKE and
${\it\epsilon}$
increase by an order of magnitude.
Dependence of the TKE (squares) and dissipation (diamonds) on particle bulk density of the two-way coupled simulations relative to a one-way coupled simulation. The increased TKE and dissipation at higher mass loading indicates the transfer of gravitational potential energy from the particles to the surrounding fluid.

This increase of turbulent agitation with particle density cannot be attributed to the loading, since the turbulent energy has been shown to decrease with increasing mass loading in stationary turbulence (Squires & Eaton Reference Squires and Eaton1991; Boivin et al.
Reference Boivin, Simonin and Squires1998). The cause is instead found in the drag forces exerted by the quickly settling particles: as a particle falls, the lost gravitational potential energy is input to the fluid. This picture is consistent with the results of the simulations of Elghobashi & Truesdell (Reference Elghobashi and Truesdell1993) for decaying turbulence as well as the experiments of Yang & Shy (Reference Yang and Shy2005) and Hwang & Eaton (Reference Hwang and Eaton2006). In particular, Elghobashi & Truesdell (Reference Elghobashi and Truesdell1993) found that the addition of heavy particles under the effect of gravity decreased the decay rate of homogeneous isotropic turbulence, and Yang & Shy (Reference Yang and Shy2005) and Hwang & Eaton (Reference Hwang and Eaton2006) found that the TKE in homogeneous turbulence in the particle-laden cases was consistent with an energy balance that takes into account the transfer of the particle potential energy to the fluid phase. A clear effect of the modification imparted to the flow structure is a strong turbulence anisotropy. This is observed in the enhancement of vertical velocity fluctuations with respect to the horizontal fluctuations, represented by the
$w_{rms}/u_{rms}$
ratio plotted in figure 4.
The ratio of the r.m.s. vertical and horizontal velocities for each of the different particle densities with two-way momentum coupling.

Figure 5 presents the energy spectra for the non-heated cases, including the one-way coupled case in which the particles do not play a role. At the relatively low density ratio of
${\it\rho}_{p}/{\it\rho}_{o}=833$
, the momentum two-way coupling injects a moderate amount of energy at the high wavenumbers and depletes energy at the low wavenumbers, which is consistent with the pivoting mechanism. The increasing density ratios increases the energy at the large scales. This is in stark contrast with previous numerical studies, in which it was found that the presence of small particles inhibits turbulence activity at low wavenumbers (e.g. Sundaram & Collins Reference Sundaram and Collins1997; Boivin et al.
Reference Boivin, Simonin and Squires1998; Bosse et al.
Reference Bosse, L. and Meiburg2006). Those studies investigated cases with either zero or moderate particle falling speed compared to the fluid velocity fluctuations (
$V_{t}/u_{{\it\eta}}\leqslant 1$
), while here the higher density ratios give
$V_{t}/u_{rms}=O(1)$
and
$V_{t}/u_{{\it\eta}}=O(10)$
. Interestingly, the intermediate scales appear to be weakly affected by the particles. The gravitational settling, which plays a major role for
${\it\rho}_{p}/{\it\rho}_{o}=4166$
and 8333, acts on the entire domain and perturbs the largest scales of fluid motion, converting the particle potential energy into turbulent energy. The dispersed phase intercepts part of the fluid energy at the intermediate scales and inject it back into the fluid at small (down to dissipative) scales. It is worth noting that Yang & Shy (Reference Yang and Shy2005), for their heavier particle case (
$St_{{\it\eta}}=1.9$
and
$V_{t}/u_{{\it\eta}}\approx 5$
), reported spectra being augmented at all scales, and more markedly at the lowest and highest wavenumbers, providing qualitative agreement with the results that we show here.
3.1.3 Settling velocity
Figure 6 shows the fractional variation of particle settling velocity with respect to the theoretical terminal velocity in quiescent air. In agreement with previous studies, the increase in settling rate ranges between 20 and 35 % (Yang & Shy Reference Yang and Shy2005; Bosse et al.
Reference Bosse, L. and Meiburg2006). The result for the one-way coupled case (
${\it\rho}_{p}/{\it\rho}_{o}=833$
,
$St_{{\it\eta}}=0.18$
) is also shown and is in excellent agreement with the results obtained in similar regimes (Wang & Maxey Reference Wang and Maxey1993; Dejoan & Monchaux Reference Dejoan and Monchaux2013; Good et al.
Reference Good, Ireland, Bewley, Bodenschatz, Collins and Warhaft2014). Depending on the normalization, the latter reads
$(W_{p}-V_{t})/u_{{\it\eta}}=0.32$
or
$(W_{p}-V_{t})/u_{rms}=0.083$
.
The turbulent energy spectra (normalized using the one-way coupled dissipation
${\it\epsilon}_{o}$
) for the different particle densities, full circles:
${\it\rho}_{p}/{\it\rho}_{o}=833$
, open circles:
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, triangles:
${\it\rho}_{p}/{\it\rho}_{o}=8333$
, crosses:
${\it\rho}_{p}/{\it\rho}_{o}=833$
with one-way coupling.

The settling parameter for the various particle densities with respect to the theoretical terminal velocity (a) and the root-mean-square velocity (b). Line: two-way coupled, cross mark: one-way coupled.

In the two-way coupled cases the increase in settling velocity is more pronounced for higher density ratios. The trend is the same whether the change in settling velocity is normalized by
$V_{t}$
or by
$u_{rms}$
. This is consistent with results of Bosse et al. (Reference Bosse, L. and Meiburg2006), who found a monotonic increase in settling rate with increasing mass loading. They described this behaviour as a consequence of the collective particle drag that enhances the downward fluid motion. However, in the present regime, another mechanism may produce this result. The heavier particles release more gravitational potential energy into the fluid, leading to substantially higher turbulence intensity (figure 3).
3.2 Settling of heated particles
3.2.1 Particle distribution
We begin our analysis of the heated cases by exploring the distribution of the particles in the turbulent flow. The PDF of Voronoi cell volumes in figure 7 confirms that the particles are strongly clustered also in this regime. Again the particles with
$St_{{\it\eta}}$
closer to unity (see table 1) display more intense clustering, as evident from the long tails of the distributions.
Voronoi volume distributions for the particles with different density ratios at different heating levels: (a)
$T_{p}/T_{o}=1.3$
, (b)
$T_{p}/T_{o}=2.2$
, full circles:
${\it\rho}_{p}/{\it\rho}_{o}=833$
, open circles:
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, triangles:
${\it\rho}_{p}/{\it\rho}_{o}=8333$
, line:
${\it\Gamma}$
distribution.

In homogeneous turbulence the mechanism of preferential concentration causes the inertial particles to be in regions of low enstrophy and high strain (Squires & Eaton Reference Squires and Eaton1991; Gibert et al.
Reference Gibert, Xu and Bodenschatz2012). This is confirmed also for the considered flow in figure 8, that displays the joint PDF of normalized particle concentration
$C/C_{o}$
and fluid enstrophy
${\it\xi}/{\it\xi}_{o}$
, where
${\it\xi}$
is the fluid enstrophy. Here
$C_{o}$
is the domain-averaged concentration and
${\it\xi}_{o}$
is the average enstrophy in the one-way coupled, unheated case. In all considered cases higher concentration levels are found in regions of low enstrophy. However, both particle density and heating have significant effects. The heavier particles (whose settling velocity largely exceeds the Kolmogorov velocity) tend to sample the flow more homogeneously. This is true despite the fact that
$Re_{{\it\lambda}}$
substantially increases with particle density (see table 1). The reason for the more uniform distribution of heavier particles may be twofold. On the one hand, the relatively high settling velocity de-correlates the particle position from the regions of high strain and low enstrophy. On the other hand, the strong momentum two-way coupling alters the structure of the turbulence, which modifies the classic picture of preferential concentration.
The logarithmically scaled joint probability density functions of normalized concentration and enstrophy. (a,d,g)
$T_{p}/T_{o}=1$
, (b,e,h)
$T_{p}/T_{o}=1.3$
, (c,f,i)
$T_{p}/T_{o}=2.2$
. (a–c)
${\it\rho}_{p}/{\it\rho}_{o}=833$
, (d–f)
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, (g–i)
${\it\rho}_{p}/{\it\rho}_{o}=8333$
. Greyscale contour levels are equally spaced (arbitrary units).
The introduction of radiative heat has a similar effect in that the particles become more uniformly distributed with respect to the enstrophy field (figure 8). Since the settling velocity is reduced by the heating (as will be shown in § 3.2.3), this is attributed to the turbulence modification. As it will be discussed in § 3.2.4, the particles release heat into the surrounding fluid, producing low density and positively buoyant plumes, which in turn modify the turbulent flow field.
In the presence of gravity, according to the mechanism of preferential sweeping proposed by Wang & Maxey (Reference Wang and Maxey1993), the inertial particles sample downward regions more often than upward regions of the turbulent flow. This is true for all considered non-heated cases, as shown in figure 9 by the joint PDF of particle concentration
$C$
and fluid vertical velocity
$w$
(normalized by
$u_{rms}$
). Also, the two-way momentum coupling causes the fluid near the heavy particles to be dragged downward even faster (Bosse et al.
Reference Bosse, L. and Meiburg2006), increasing the probability of finding particles in regions
$w<0$
. However, in the presence of heating, the trend is reversed, and the particles are more likely to be found in regions of upward turbulent fluctuations. This is again a hint of the presence of positively buoyant plumes shed by the hot particles. It is also apparent that for the non-heated cases the range of vertical fluid velocities
$w$
becomes broader for more massive particles, which reflects the increased turbulence activity (see figure 3). Instead, with increasing radiation, the range of
$w$
does not broaden monotonically, which is consistent with the trend of TKE to be discussed in § 3.2.4.
The logarithmically scaled joint probability density functions of normalized concentration and vertical velocity: (a–c)
$T_{p}/T_{o}=1$
, (d–f)
$T_{p}/T_{o}=1.3$
, (g–i)
$T_{p}/T_{o}=2.2$
. (a,d,g)
${\it\rho}_{p}/{\it\rho}_{o}=833$
, (b,e,h)
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, (c,f,i)
${\it\rho}_{p}/{\it\rho}_{o}=8333$
. Greyscale contour levels are equally spaced (arbitrary units).
3.2.2 Statistical evidence of buoyant plumes
Given that particles are the sole heat source for the gas, we expect positive correlation between gas temperature and particle density. This is confirmed in figure 10, showing the joint PDF of particle concentration and local fluid density.
The logarithmically scaled joint probability density functions of normalized concentration and fluid density: (a,c,e)
$T_{p}/T_{o}=1.3$
, (b,d,f)
$T_{p}/T_{o}=2.2$
. (a,b)
${\it\rho}_{p}/{\it\rho}_{o}=833$
, (c,d)
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, (e,f)
${\it\rho}_{p}/{\it\rho}_{o}=8333$
. Greyscale contour levels are equally spaced (arbitrary units).
While heating the particles broadens the range of fluid density, the joint PDF of
$w/u_{rms}$
and
${\it\rho}/{\it\rho}_{o}$
in figure 11 indicates that a correlation persists between the low density regions and the upward velocity fluctuations. This suggests that, even in the highly turbulent environment through which the heavy particles settle, hot regions of upward fluid motion remain a distinct feature.
The joint probability density functions of normalized fluid velocity and fluid density: (a,c,e)
$T_{p}/T_{o}=1.3$
, (b,d,f)
$T_{p}/T_{o}=2.2$
. (a,b)
${\it\rho}_{p}/{\it\rho}_{o}=833$
, (c,d)
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, (e,f)
${\it\rho}_{p}/{\it\rho}_{o}=8333$
. Greyscale contour levels are equally spaced (arbitrary units).

3.2.3 Settling velocity
Figure 12 presents the difference between the average particle settling velocity and the nominal terminal velocity as a function of the particle-to-fluid temperature ratio. In order to highlight the impact on the particle–turbulence interaction, we normalize the velocity difference by
$u_{rms}$
(normalization using
$u_{{\it\eta}}$
as the fluid velocity scale leads to a similar result). For all density ratios the heating of the particles hinders the falling speed. We remark that this is not a consequence of a global upward fluid velocity, since the mean velocity is set to zero at every time step. It is rather caused by the high temperature regions around the particles (and especially the particle clusters), which produce localized updraft in correspondence of the high particle concentration areas. At the most intense radiation level this effect offsets the preferential sweeping mechanism, and the falling speed is lower than the nominal terminal velocity in a quiescent fluid.
Plot of the settling parameter for each of the different particle densities for different heating levels and different density ratios, full circles:
${\it\rho}_{p}/{\it\rho}_{o}=833$
, open circles:
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, triangles:
${\it\rho}_{p}/{\it\rho}_{o}=8333$
.

Interestingly, the heavier particles reduce their falling speed more than the lighter particles. A possible reason is that the light particles are less densely clustered (figure 7). However, this cannot be the sole explanation for the observed behaviour, since the particles with
${\it\rho}_{p}/{\it\rho}_{o}=4166$
are slightly more clustered than the particles with
${\it\rho}_{p}/{\it\rho}_{o}=8333$
. An additional factor is that the buoyancy production term in the turbulent kinetic energy transport is larger for particles of greater mass (as discussed in § 3.2.4 below). This is associated with more extended plumes of positively buoyant fluid around the particle clusters, and therefore a stronger reduction of settling velocity.
The PDF of the particle settling velocity for each case is shown in figure 13. The PDFs are centred about the nominal settling velocity of the particle in a quiescent fluid and are normalized. Thus if the mean is negative, the particles tend to experience settling enhancement due to turbulence. If the mean is positive, the particles fall slower than in a quiescent fluid. All distributions have a nearly Gaussian shape, with the position of the mean reflecting the trend reported in figure 12: heavier particles experience a stronger hindering of the settling when heated. The change of the variance with particle heating is not intuitive. For
${\it\rho}_{p}/{\it\rho}_{o}=833$
, the PDF broadens with radiation intensity, while for heavier particles the PDF is first narrowed and then broadened by increasing radiation. This non-monotonic behaviour, to which we alluded in § 3.2.1, is the consequence of the non-monotonic change in TKE with radiation intensity, which is discussed below (§ 3.2.4).
The probability density functions of the particle settling velocity normalized by the nominal terminal velocity for
${\it\rho}_{p}/{\it\rho}_{o}=833$
(a),
${\it\rho}_{p}/{\it\rho}_{o}=4166$
(b), and
${\it\rho}_{p}/{\it\rho}_{o}=8333$
(c). Plus signs:
$T_{p}/T_{o}=1$
, squares:
$T_{p}/T_{o}=1.3$
, diamonds:
$T_{p}/T_{o}=2.2$
.

3.2.4 Energy transfer to the fluid phase
As already demonstrated, the presence of inertial particles has a substantial impact on the turbulence statistics. In the following sections we further highlight the turbulence modification when the particles are heated.
To assess the contribution of the heating to the flow energy budget, we consider the transport equation for the turbulent kinetic energy
$$\begin{eqnarray}\displaystyle \frac{\partial k}{\partial t}+u_{j}\frac{\partial k}{\partial x_{j}} & = & \displaystyle \frac{1}{{\it\rho}}\frac{\partial (u_{i}{\it\tau}_{ij})}{\partial x_{j}}+{\it\rho}g_{i}u_{i}-{\it\epsilon}-\frac{1}{{\it\rho}}\frac{\partial (pu_{i})}{\partial x_{i}}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1}{{\it\rho}}\frac{\partial u_{i}}{\partial x_{i}}p-\frac{1}{{\it\rho}}\mathop{\sum }_{n=1}^{N_{p}}\frac{m_{p}}{{\it\tau}_{p}}(2k-u_{i}v_{i,n}){\it\delta}(\boldsymbol{x}-\boldsymbol{y}_{n}),\end{eqnarray}$$
where the dissipation
${\it\epsilon}$
for a Newtonian fluid with zero bulk viscosity is
and, as a counterpart, the equation for the kinetic energy of a Lagrangian particle
$k_{p}$
:
The terms in the right-hand side of the turbulent kinetic energy are the transport terms (including viscous diffusion and pressure strain), the viscous dissipation, and the particle–fluid interaction. The particle kinetic energy equation has an opposite particle–fluid interaction term and a gravity production term. The ratio of the buoyancy production (
${\it\Pi}_{b}=-({\it\rho}-{\it\rho}_{o})gw$
) to the dissipation (
${\it\epsilon}$
) is plotted in figure 14 as a function of the particle-to-fluid temperature ratio. As the heating level increases, the local density fluctuations become more intense and the buoyancy forces account for a larger fraction of the turbulent kinetic energy production. The effect of the particle density is also apparent: denser particles possess (under the present assumption of constant material properties) higher thermal capacity and therefore can store more energy and sustain higher temperatures against convection losses.
Ratio of the buoyancy production to the viscous dissipation for the different particle densities as a function of heating level, full circle:
${\it\rho}_{p}/{\it\rho}_{o}=833$
, open circle:
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, triangle:
${\it\rho}_{p}/{\it\rho}_{o}=8333$
.

The balance of the various terms in TKE transport is depicted in figure 15, where the ratio of
$\text{TKE}/\text{TKE}_{{\it\Phi}=0}$
is plotted as a function of the particle temperature for the different density ratios. When looking at the effect of particle heating, a non-monotonic trend is apparent: TKE first decreases at intermediate particle temperatures, and then increases at higher temperatures. Such a trend would not be visible for the lighter particles (
${\it\rho}_{p}/{\it\rho}_{o}=833$
) if we had not added case at
$T_{p}/T_{o}=1.006$
. The latter shows a slight decrease in TKE before increasing at higher temperature ratios. This suggests that a minimum exists for all radiation levels, and that the temperature ratio that minimizes TKE is higher for heavier particles. A denser matrix of cases is needed to quantitatively characterize this behaviour, but the present set of simulations is sufficient to demonstrate this qualitatively.
The TKE for each particle density as a function of heating level, including the low heating case for
${\it\rho}_{p}/{\it\rho}_{0}=833$
, full circles:
${\it\rho}_{p}/{\it\rho}_{o}=833$
, open circles:
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, triangles:
${\it\rho}_{p}/{\it\rho}_{o}=8333$
.

The origin of this non-monotonic dependence of TKE on the particle heating is to be found in the competing effects of settling and buoyancy, which inject momentum into the flow in opposite directions. The settling heavy particles transfers potential energy to the fluid through drag force. On the other hand, the radiation-induced buoyancy reduces the settling velocity (and more so for heavier particles) which indirectly reduces TKE. However the buoyancy term is a source of turbulence production, and at sufficiently high heating levels its positive contribution to the turbulence dominates, causing TKE to grow.
The impact of particle heating on the fluid motion at the various scales is illustrated by the energy spectra in figure 16. It appears that the radiation, both when reducing and when increasing the turbulence activity, affects the intermediate and especially the large scales of fluid motion. This is consistent with the idea that the change in TKE results from the balance of settling velocity and buoyancy production. Indeed, gravitational settling affects mostly the large scales, while buoyancy-induced updrafts occur at the scales of the particle clusters, which are expected to have sizes ranging between the intermediate and large scales of the turbulence. It is also noteworthy that the magnitude of the spectra show the same non-monotonicity with respect to radiative input that the TKE does.
The ratio of the two-way coupled energy spectra to the one-way coupled energy spectrum for each of the heating levels and particle densities as a function of the normalized wavenumber. (a)
${\it\rho}_{p}/{\it\rho}_{o}=833$
, (b)
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, (c)
${\it\rho}_{p}/{\it\rho}_{o}=8333$
. Plus signs:
$T_{p}/T_{o}=1$
, squares:
$T_{p}/T_{o}=1.3$
, diamonds:
$T_{p}/T_{o}=2.2$
, crosses: one-way coupled,
$T_{p}/T_{o}=1$
.

3.2.5 Turbulence anisotropy
The downward drag force exerted by the falling particles on the fluid modifies the turbulent motion at all scales, as suggested by the spectra in figure 16. The force also breaks the isotropy of the turbulence, as demonstrated by figure 17, which plots the longitudinal integral scales (determined from the corresponding two-point correlations) in the vertical and horizontal directions,
$L_{w,zz}$
and
$L_{u,xx}$
, respectively, defined by the distance at which the velocity two-point correlation first reaches zero. The presence of the particles increases the integral length scale in the direction of gravity, and more so for the heavier particles. The effect of heating makes this effect somewhat more pronounced, but the qualitative trend is unchanged. On the other hand, in the horizontal direction, the addition of particles reduces the integral length scale in the absence of thermal effects. This is consistent with the redistribution of turbulent energy from the horizontal to the vertical direction (Elghobashi & Truesdell Reference Elghobashi and Truesdell1993). However, when heating the particles, the horizontal length scale increases, indicating the level of anisotropy is somewhat reduced. The effect is confirmed by figure 18, which shows the ratio of the vertical to horizontal components of the velocity fluctuations as a function of particle heating. The trend demonstrates that particle heating augments the anisotropy in the
${\it\rho}_{p}/{\it\rho}_{o}=833$
case, but reduces it for the higher density ratios. This is because, as discussed in § 3.2.3, heating reduces the falling speed of heavier particles more sharply. This limits the particle momentum transfer to the flow, and in turn lessens anisotropy and turbulence enhancement. Thus the particle heating acts to reduce the settling speed and redistributes the turbulent energy from the vertical to the horizontal directions.
The vertical and horizontal velocity two-point correlation length scales measured with respect to the box size
$L_{o}$
as a function of particle density for each of the heating levels, plus signs:
$T_{p}/T_{o}=1$
, squares:
$T_{p}/T_{o}=1.3$
, diamonds:
$T_{p}/T_{o}=2.2$
.

Ratio of the root-mean-square vertical and horizontal velocities as a function heating level for each of the different particle densities, full circles:
${\it\rho}_{p}/{\it\rho}_{o}=833$
, open circles:
${\it\rho}_{p}/{\it\rho}_{o}=4166$
, triangles:
${\it\rho}_{p}/{\it\rho}_{o}=8333$
.

4 Conclusions
We investigated solid particles falling through homogeneous air turbulence while heated by thermal radiation. In the considered optically thin regime, radiation is simply modelled as a fixed heat flux input to each particle. Particles much smaller than the Kolmogorov scale are considered, which permits the solution for all relevant flow scales taking into account both the momentum and thermal coupling between phases. In the considered range of particle density and heating levels, the Stokes number indicates strong preferential concentration and the settling velocity is comparable to the fluid velocity fluctuations. The considered mass loading is sufficient to cause sizable modification of the air turbulence, and the maximum radiation level more than doubles the absolute temperature of the particles. The main findings can be summarized as follows.
-
(1) Clustering is intense for particles with
$St_{{\it\eta}}=O(1)$
, regardless of momentum or temperature coupling. However, the mechanism of preferential concentration in high-strain/low-vorticity regions is blurred by the turbulence modification. -
(2) In the absence of radiation, particles with settling parameter
$V_{t}/u_{{\it\eta}}\approx 1$
experience substantial increase in their settling velocity due to preferential sweeping, as expected from previous studies. Settling of heavier particles (
$V_{t}/u_{{\it\eta}}>1$
) augments the turbulent energy and increases the vertical fluid velocity fluctuations with respect to the horizontal component. -
(3) In the presence of radiation, hot buoyant plumes are shed from the heated particles, which reduces the settling velocity. The preferential sweeping mechanism breaks down and particles do not preferentially sample regions of downward-moving fluid. When normalized with respect to the turbulent velocity fluctuations, the settling velocity is more strongly reduced for heavier particles.
-
(4) Particle heating reduces turbulent kinetic energy at low radiation levels, and increases it at high radiation levels, with the minimum depending on the particle density. This non-monotonic trend is the consequence of two competing effects of the radiative forcing: reducing the settling velocity (which reduces turbulence) and causing buoyancy-driven motion (which enhances turbulence).
This is the first study addressing the settling of heated solid particles in turbulence. Both the heated and the non-heated regimes we investigated occupy positions in the parameter space which had not been explored in the past, despite being highly relevant for a wealth of natural and environmental problems. Although the design of systems containing heated particles is often assumed to have decoupled thermal and momentum effects, we have shown here that the coupling between the two processes is very strong and has a non-trivial effect on the flow statistics. While we identified the mechanisms of particle–turbulence interaction, further research is warranted to better characterize the individual role of the various parameters at play, including mass loading, settling velocity, Stokes number and temperature ratio, as well as to extend the physics to include more complex phenomena such as evaporation, chemistry and radiative scattering and blackbody emission.
Acknowledgements
This work was supported by the United States Department of Energy under the Predictive Science Academic Alliance Program 2 (PSAAP2) at Stanford University. We are grateful for the insightful discussions with S. Subramaniam and H. Kolla for bringing their insight during the CTR Summer Program.
Appendix A. Non-dimensional equations
To illustrate the roles of various terms in the governing equations, we present here the non-dimensionalized governing equations. We define the non-dimensional variables as
where the subscript
$o$
denotes a reference value. Under these scalings, the Navier–Stokes equations become
$$\begin{eqnarray}\displaystyle \frac{\partial ({\it\rho}^{\ast }u_{i}^{\ast })}{\partial t^{\ast }}+\frac{\partial ({\it\rho}^{\ast }u_{i}^{\ast }u_{j}^{\ast })}{\partial x_{j}^{\ast }} & = & \displaystyle -\frac{\partial p^{\ast }}{\partial x_{i}^{\ast }}+\frac{1}{Re}\left(\frac{\partial ^{2}u_{i}^{\ast }}{\partial x_{j}^{\ast }\partial x_{j}^{\ast }}+\frac{1}{3}\frac{\partial }{\partial x_{i}^{\ast }}\frac{\partial u_{k}^{\ast }}{\partial x_{k}^{\ast }}\right)-\frac{1}{Fr^{2}}{\it\rho}^{\ast }{\it\delta}_{i3}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{{\it\phi}_{m}}{N_{p}St}\mathop{\sum }_{n=1}^{N_{p}}(u_{i}^{\ast }-v_{i,n}^{\ast }){\it\delta}(\boldsymbol{x}^{\ast }-\boldsymbol{y}_{n}^{\ast }),\end{eqnarray}$$
where
$Re={\it\rho}_{o}u_{o}L_{o}/{\it\mu}$
is the Reynolds number,
$Fr=u_{o}/\sqrt{gL_{o}}$
is the Froude number,
${\it\phi}_{m}=m_{p}N_{p}/{\it\rho}_{o}L_{o}^{3}$
is the mass loading ratio for a particle mass
$m_{p}$
,
$N_{p}$
is the number of particles and
$St$
is the Stokes number. Note that we have substituted in the definition of the viscous stress tensor
${\it\tau}_{ij}={\it\mu}((\partial u_{i}/\partial x_{j})+(\partial u_{j}/\partial x_{i})-(2/3)(\partial u_{k}/\partial x_{k}){\it\delta}_{ij})$
. The energy equation for the fluid phase is given by
where
$Pr={\it\mu}C_{p}/k$
is the Prandtl number,
${\it\gamma}=C_{p}/C_{v}$
is the adiabatic constant and
$Nu$
is the Nusselt number of the particles. Under the assumption of an ideal gas with a low-Mach flow, the ideal gas law reduces to
The non-dimensionalized Lagrangian particle equations are
where
${\it\Gamma}=C_{v,p}/C_{v}$
is the ratio of the particle specific heat capacity to the fluid isochoric specific heat capacity, and
$I_{o}^{\ast }=I_{o}D/4kT_{o}$
is the non-dimensional radiation source.
In this study we varied the particle density
${\it\rho}_{p}$
, and the radiation intensity
$I_{o}$
. The non-dimensional radiation source and mass loading ratio are thus directly varied. The Reynolds, Stokes and Froude numbers all have a non-trivial dependence on the velocity and length scales in the turbulence and hence, due to the momentum and thermal coupling, all change in response to the variation in the particle density and radiation intensity. In all cases, the Nusselt number is set equal to
$2$
by the assumption of Stokesian convection from each particle. The Prandtl number and adiabatic constants are also assumed to be constant for the fluid and set equal to values typical for air at atmospheric conditions. The number of particles in each case was also held constant.
















