1 Introduction
Particle-laden turbulent flows are present in a broad range of natural and industrial phenomena, including the formation of clouds, urban pollutant dispersion, protoplanetary disks, volcanic eruptions and solid fuel combustion. In the simplest case, the particles are passive, and are purely advected by the flow. Tracer-based measurement techniques such as particle image velocimetry (PIV) rely on such an assumption. However, in many cases there is a complex two-way coupling where mass, momentum and energy of both the particles and the fluid phase are affected by each other. This is even more complicated if the fluid phase is turbulent since the particles can fundamentally disrupt the multi-scale turbulence cascade. These complexities make the study of particle-laden turbulent flows challenging.
The specific interests of this study are scenarios with significant heat transfer between the dispersed phase and fluid phase. The heat transfer can be sustained either by radiative heating of the particles (Zamansky et al. Reference Zamansky, Coletti, Massot and Mani2014) or by an internal exothermic reaction in particles (Mashayek Reference Mashayek2000). Therefore, in this study we assume particles are heated through an external source (e.g. radiation), and then conductively transfer heat to the gas phase.
Even in the absence of particle heating, there is a range of well-studied phenomena that are caused by the particle interaction with the turbulent eddies, such as preferential concentration of particles in regions of low vorticity (Squires & Eaton Reference Squires and Eaton1991), preferential sweeping of particles with high inertia in cases where gravity is present (Wang & Maxey Reference Wang and Maxey1993) and fluid acceleration by particle drag at high particle concentrations which in turn enhances particle settling velocity (Bosse, Kleiser & Meiburg Reference Bosse, Kleiser and Meiburg2006). In figure 1 the particle distribution in particle-laden homogeneous isotropic turbulence is shown for various Stokes numbers (the ratio of particle relaxation time to the Kolmogorov time). Evidently, for Stokes number close to unity, particles are preferentially concentrated. In §§ 4.2 and 5, using radial distribution functions, we quantify the preferential concentration of particles.
Figure 1. Particle distribution in a particle-laden homogeneous isotropic turbulence simulation. Only particles within a tiny slice are presented. The Stokes number is based on Kolmogorov time scale.
External heating of preferentially concentrated particles results in inhomogeneous conductive heat transfer from the dispersed phase to the gas phase. The gas phase, consequently, experiences inhomogeneous expansion at the location of particle clusters. In this work, we study the effects of preferential expansion of the gas field on the spectral properties of the turbulent flow. We study the effect of particle clustering in § 5 by considering systems over a wide range of Stokes numbers.
A relevant application of this study is the particle-based solar receiver (Pouransari & Mani Reference Pouransari and Mani2017), in which solid particles dispersed in background fluid are thermally heated by solar radiation. The fluid is radiatively transparent; however, it indirectly absorbs radiation energy via conductive heat transfer from hot solid particles. The downstream hot fluid can then be employed in some power generation cycle, e.g. Brayton cycle (Mancini, Kolb & Prairie Reference Mancini, Kolb and Prairie1997), or in methane splitting reactors (Kogan & Kogan Reference Kogan and Kogan2003). This idea was first proposed by Abdelrahman, Fumeaux & Suter (Reference Abdelrahman, Fumeaux and Suter1979) and Hunt (Reference Hunt1979). The optimal design of such receivers requires a comprehensive understanding of the coupling between the heating of the particles and the dynamics of the carrier fluid. In particular, we investigate the effects of particle heating on the dynamics of the carrier turbulent flow.
The results presented in this work may also be found relevant in other systems involving heat exchange between a dispersed phase and a gas, such as volcanic eruptions, circumstellar disks, injection of fuel spray in a combustion chamber and rocket engines. Conceivably, heated particles introduce additional coupling phenomena such as fluid dilatation and velocity fluctuations due to particle-to-gas heat transfer. In an analogous setting (turbulent combustion), flame-induced fluid dilatation was observed to affect the turbulent kinetic energy (TKE) spectrum at wavenumbers corresponding to the flame scales (Kolla et al. Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014). We are interested in investigating whether a similar phenomenon arises in heated particle-laden flows. In this work, numerical simulations are used to study the kinetic energy spectra of the homogeneous isotropic turbulent (HIT) flows with heated particles. The dispersed phase is modelled as Stokesian point particles. We compare the energy spectra of a cold HIT with a heated particle-laden HIT. Unlike the turbulent combustion setting considered by Kolla et al. (Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014), in our setting all three spatial directions are statistically homogeneous. This property allows rigorous investigation of turbulence in terms of wavenumbers and its decomposition to dilatational and divergence-free modes.
Utilizing our density-weighted formulation, we study flow dilatation and its influence on the energy cascade across various scales. The energy transfer analysis reveals any presence of backscatter of the flow kinetic energy. Essentially, despite the conventional understanding about the HIT where energy cascades from low to high wavenumbers, here energy is injected into the fluid at wavenumbers corresponding to the length scales of particle clusters through flow dilatation. This is an important consideration in any subgrid-scale (SGS) modelling and large-eddy simulation (LES) of turbulent flows involving inhomogeneous heating.
2 Mathematical framework
A novel aspect of this study is the use of a variable-density formulation to study kinetic energy spectra. Variations in the density field modify turbulence at scales comparable to the Kolmogorov scales. This is the case in both turbulent combustion and the present application. The traditional mathematical framework for studying spectra and energy dynamics in spectral space is valid for constant-density flows and cannot be applied in a straightforward fashion to variable-density flows. Kolla et al. (Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014) proposed a density-weighted formulation for two-point velocity correlation tensors that allows the study of energy spectra in variable-density flows. The formulation, based on conservation laws for variable-density flows, yields a balance equation for the Favre-averaged TKE in wavenumber space. The balance equation illuminates the role of various quantities such as pressure fluctuation, dilatation, mean velocity gradient, density fluctuation and dissipation on the kinetic energy balance in wavenumber space.
The framework is centred on the following definition for the two-point velocity correlation tensor,
$$\begin{eqnarray}\widetilde{\unicode[STIX]{x1D619}}_{ij}(\boldsymbol{r})={\textstyle \frac{1}{2}}\langle \unicode[STIX]{x1D70C}u_{i}(\boldsymbol{x})u_{j}(\boldsymbol{x}+\boldsymbol{r})+u_{i}(\boldsymbol{x})\unicode[STIX]{x1D70C}u_{j}(\boldsymbol{x}+\boldsymbol{r})\rangle ,\end{eqnarray}$$
where
$\unicode[STIX]{x1D70C}$
is density,
$u_{i}$
denotes the
$i$
th component of velocity fluctuation,
$\langle ~\rangle$
denotes an averaging over the position vector
$\boldsymbol{x}$
in homogeneous directions, and the vector
$\boldsymbol{r}$
denotes the separation vector. Note that, upon contracting the tensorial indices and setting the separation vector to zero, the two-point correlation tensor yields the Favre-averaged TKE, analogous to its constant-density counterpart. Using this property and starting from conservation equations for
$\unicode[STIX]{x1D70C}u_{i}$
and
$u_{i}$
in variable-density flows, a balance equation for the Favre-averaged TKE in wavenumber space was derived and presented in Kolla et al. (Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014). For the homogeneous isotropic configuration chosen in the present study, there are no gradients of mean quantities, and thus only three terms remain non-zero in the energy spectral balance equation (Kolla et al.
Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014): a pressure-dilatation term that represents a source, a density-weighted triple velocity correlation term that transfers energy across wavenumbers, and a molecular dissipation term. Here, we focus on the first two terms.
The pressure-dilatation term encapsulates the effect of fluid dilatation due to particle heating. This term is of the form
$$\begin{eqnarray}-\left\langle \text{Re}\left\{{\mathcal{F}}^{\ast }\{u_{j}\}{\mathcal{F}}\left\{\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{j}}\right\}\right\}+\text{Re}\left\{{\mathcal{F}}^{\ast }\{\unicode[STIX]{x1D70C}u_{j}\}{\mathcal{F}}\left\{\frac{1}{\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{j}}\right\}\right\}\right\rangle ,\end{eqnarray}$$
where
${\mathcal{F}}$
,
${\mathcal{F}}^{\ast }$
and
$\text{Re}$
denote the Fourier transform, the complex conjugate of the Fourier transform and the real part, respectively. The pressure-dilatation term denotes the primary contribution of density variation through dilatation in the TKE.
The density-weighted triple velocity correlation term represents energy transfer among wavenumbers. Note that the integral of the energy transfer term over all wavenumbers is identically zero. This term is of the form
$$\begin{eqnarray}\left\langle \text{Re}\left\{{\mathcal{F}}^{\ast }\{\unicode[STIX]{x1D70C}u_{j}\}{\mathcal{F}}\left\{u_{k}\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{k}}\right\}\right\}+\text{Re}\left\{{\mathcal{F}}^{\ast }\{u_{j}\}{\mathcal{F}}\left\{\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}u_{j}u_{k}}{\unicode[STIX]{x2202}x_{k}}\right\}\right\}\right\rangle .\end{eqnarray}$$
The significance of the energy transfer term is that it quantifies the degree of backscatter that might arise from the energy injected at high wavenumbers by the pressure-dilatation term. In (2.2) and (2.3),
$\langle ~\rangle$
denotes an averaging in the wave space.
3 Simulation framework
3.1 Governing equations
The flow conditions correspond to the low-Mach-number regime with a constant thermodynamic pressure in space, which is governed by the ideal gas law. In this study the nominal Mach number, defined based on turbulent fluctuation velocity, is
$1.8\times 10^{-3}$
, which is much smaller than unity. Therefore, we can safely neglect the contribution of local hydrodynamic pressure in the thermodynamic pressure that is used in the equation of state. Variation in the gas density field results in a variable-coefficient Poisson equation for the hydrodynamic pressure. The Poisson equation is solved using an iterative spectral method. The details of the numerical methods are explained in Pouransari, Mortazavi & Mani (Reference Pouransari, Mortazavi and Mani2015). Zamansky et al. (Reference Zamansky, Coletti, Massot and Mani2014) showed that a particle-laden turbulence subjected to external radiation and gravity produces a sustained turbulence. We ignore the effect of gravity by considering small domains and fast processes (large Froude number). This is justifiable in regimes where particles are heated in a stream, such as particle-based solar receivers, where the residence time is much faster than the time scale of buoyancy. In addition, with this simplification we can compare the effect of dilatation in the heated particle-laden problem and analogous regimes in combustion.
With the above assumptions, the gas phase can be described with the following set of non-dimensionalized conservation equations:
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(\unicode[STIX]{x1D70C}u_{j})=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}u_{i})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(\unicode[STIX]{x1D70C}u_{i}u_{j})=-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x_{i}}+\frac{1}{\mathit{Re}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}-\frac{2}{3}\frac{\unicode[STIX]{x2202}u_{k}}{\unicode[STIX]{x2202}x_{k}}\unicode[STIX]{x1D6FF}_{ij}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{\unicode[STIX]{x1D6FE}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70C}T_{g})+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x_{j}}(\unicode[STIX]{x1D70C}T_{g}u_{j})=\frac{1}{\mathit{Re}\mathit{Pr}}\frac{\unicode[STIX]{x2202}^{2}T_{g}}{\unicode[STIX]{x2202}x_{j}\unicode[STIX]{x2202}x_{j}}+Q. & \displaystyle\end{eqnarray}$$
In (3.3),
$Q$
is the non-dimensional rate of heat transfer from the particles to the gas phase:
$$\begin{eqnarray}Q=\frac{D_{p}}{L}\frac{\unicode[STIX]{x03C0}\mathit{Nu}}{\mathit{Re}\mathit{Pr}}\mathop{\sum }_{i=1}^{N_{p}}(T_{p_{i}}-T_{g})\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{p_{i}}),\end{eqnarray}$$
where
$D_{p}$
is particle diameter,
$L$
is the domain size,
$T_{p_{i}}$
is the temperature of the
$i$
th particle located at
$\boldsymbol{x}_{p_{i}}$
, and
$\unicode[STIX]{x1D6FF}$
is the non-dimensional Dirac delta function. The summation term in the above equation is numerically approximated using trilinear projection/interpolation with normalization condition that the discrete integral of
$\unicode[STIX]{x1D6FF}$
over the domain is 1.
The turbulence is initialized with
$\mathit{Re}=Lu_{rms}/\unicode[STIX]{x1D708}=1550$
, where
$u_{rms}$
is the root mean square of the initial velocity field in one direction, and
$\unicode[STIX]{x1D708}$
is the kinematic viscosity. This condition corresponds to the Reynolds number based on the Taylor microscale
$\mathit{Re}_{\unicode[STIX]{x1D706}}=1.7\mathit{Re}^{1/2}=67$
(see e.g. Rosales & Meneveau Reference Rosales and Meneveau2005). Considering a typical gas-phase medium, we assume that the Prandtl number and heat capacity ratio are
$\mathit{Pr}=0.7$
and
$\unicode[STIX]{x1D6FE}=1.4$
, respectively.
In this study, we choose a set of parameters for particles relevant to particle-based solar receivers similar to those of Frankel et al. (Reference Frankel, Pouransari, Coletti and Mani2016). Considering small particles with negligible thermal Péclet number, we set
$\mathit{Nu}=2$
. The nominal (base) number of particles is
$N_{p}=1.28\times 10^{6}$
. Additionally, we present systems in which the number of particles is varied. In the nominal case,
$d_{p}/L=3\times 10^{-4}$
. This corresponds to a highly dilute case with small volume fraction of
$\unicode[STIX]{x1D6F7}=1.8\times 10^{-5}$
, and small particle size compared to the Kolmogorov length scale.
The ratio of particle relaxation time to the Kolmogorov time (Stokes number) is
$St_{\unicode[STIX]{x1D702}}=5.3$
at the initial time. As turbulence decays, the Stokes number decreases to
$O(1)$
(see figure 8). We intentionally considered this Stokes number to have high particle preferential concentration. This provides us with strong inhomogeneous heat transfer from the dispersed phase to the gas phase. In § 5, we study additional cases with different Stokes number.
Note that the ratio of particle thermal relaxation time to momentum relaxation time is
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D70F}_{th}}{\unicode[STIX]{x1D70F}_{p}}=\frac{3\mathit{Pr}\unicode[STIX]{x1D6E4}}{\mathit{Nu}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D6E4}$
is the ratio of particle specific heat capacity to that of a gas at constant pressure. In the nominal scenario, we set this ratio to 0.42, considering a typical gas–solid system. Therefore, particle-to-gas heat transfer occurs fast enough before the turbulence decays.
The energy of a particle is governed by the following non-dimensional equation:
$$\begin{eqnarray}\unicode[STIX]{x1D70F}_{th}^{\ast }\frac{\text{d}T_{p}}{\text{d}t}={\mathcal{R}}-T_{p}+{\mathcal{I}}(T_{g}),\end{eqnarray}$$
where
${\mathcal{I}}$
is the interpolation operator. We used a trilinear interpolation to find the temperature of the gas at the location of the particle. In (3.6),
$\unicode[STIX]{x1D70F}_{th}^{\ast }$
is the dimensionless particle thermal relaxation time, which can be obtained using (3.5) and the Stokes number.
Also in (3.6),
${\mathcal{R}}=R/(\mathit{Nu}\unicode[STIX]{x03C0}D_{p}kT_{0})$
is the ratio of the external heating rate received by each particle,
$R$
, to the heat transfer rate from a particle to the gas per unit temperature difference. The nominal value of
${\mathcal{R}}$
is
$0.34$
. In order to keep the effective
${\mathcal{R}}$
constant in the window of our calculation, we assume that thermodynamic pressure is constant in time. This is equivalent to a system with heating proportional to the mean temperature of the system. We also present additional calculations in which we vary this parameter systematically.
The described dimensionless parameters provide a complete set that fully determines the system under investigation. One important dimensionless group, which was not explicit in the governing equations, is
$\unicode[STIX]{x1D6FC}=\unicode[STIX]{x1D6FE}N_{p}{\mathcal{R}}\times (D_{p}/L)\times (\unicode[STIX]{x03C0}\mathit{Nu}/\mathit{Re}\mathit{Pr})$
. In terms of dimensional input parameters, this expression can be written as
$\unicode[STIX]{x1D6FC}=N_{p}R/(\unicode[STIX]{x1D70C}_{0}C_{v}L^{3}T_{0})\times (L/u_{rms})$
, which is the ratio of eddy time scale
$O(L/u_{rms})$
to the time scale of increase in internal energy due to heating. For the nominal condition considered here, this ratio is 1.0. Cases that involve heating rates much slower than the large-eddy time scale are not expected to involve significant turbulence modification due to heating. Here we investigate scenarios where the heating time scale is equal to or faster than the large eddy time.
Given the small volume fraction, the considered regimes lie in the category of dilute suspension (Balachandar & Eaton Reference Balachandar and Eaton2010). The only significant term in the dynamics of the particle equation is the Stokes drag term (Gatignol Reference Gatignol1983; Maxey & Riley Reference Maxey and Riley1983). For simplicity, the thermodynamic properties are considered temperature-independent. Gas and particles are one-way coupled for momentum and two-way coupled for energy.
Figure 2. The averaged normalized TKE as a function of time in the forced HIT. Red plots show evolution of the normalized TKE for decaying turbulence and active particle heating.
3.2 Simulation scenarios
Numerical simulation of a decaying HIT laden with particles is performed using a Lagrangian framework for particle tracking, and an Eulerian framework for simulation of the gas phase in a triply periodic domain resolved on a uniform staggered mesh similar to Mashayek & Jaberi (Reference Mashayek and Jaberi1999). Our numerical scheme is second-order-accurate in space, and fourth-order-accurate in time for both gas and dispersed phases.
Each simulation is initialized with a random turbulent flow field similar to Kwak, Reynolds & Ferziger (Reference Kwak, Reynolds and Ferziger1975). Cold particles are initiated randomly with uniform distribution in the domain. Turbulence is sustained using artificial linear forcing described in Rosales & Meneveau (Reference Rosales and Meneveau2005). After initialization, the simulation is continued without external heating to the particles until the turbulence and spatial distribution of particles become fully developed. When the particle-laden turbulence is fully developed, we identify snapshots in time for which the spatially averaged TKE,
$\langle \unicode[STIX]{x1D70C}u^{2}\rangle$
, is equal to its long-term averaged value,
$\langle \langle \unicode[STIX]{x1D70C}u^{2}\rangle \rangle$
. Note that, owing to the homogeneity,
$u$
can be any component of the velocity field. In figure 2 the evolution of the normalized TKE as a function of large-eddy turnover time is plotted during the forcing period. The large-eddy turnover time is defined as
$u_{rms}^{2}/\unicode[STIX]{x1D716}$
, where
$\unicode[STIX]{x1D716}$
is the initial dissipation rate. The circles correspond to instances when the TKE of the instantaneous field is equal to the long-term averaged TKE. Each one of these snapshots is used as the initial condition for a target calculation in which turbulence decays (forcing turned off) and particle heating is present. For improved convergence, multiple realizations are used for each case to obtain averaged statistics. Each red curve in figure 2 corresponds to one such realization. Hence, the particle heating intensity introduced in (3.6) jumps from 0 to a finite value
${\mathcal{R}}$
at the beginning of each realization (i.e. red curves in figure 2).
4 Results
We present the results that are obtained by ensemble averaging over five independent realizations. The one-dimensional TKE spectrum for a variable-density flow is obtained as follows:
$$\begin{eqnarray}E_{k}=\iint \text{Re}\{{\mathcal{F}}^{\ast }\{\unicode[STIX]{x1D70C}u_{j}\}{\mathcal{F}}\{u_{j}\}\}.\end{eqnarray}$$
In the above equation the double integral is taken over a spherical shell in the wave space.
As mentioned in § 3, the gas-phase solution points are stored on a uniform staggered grid, i.e. momentum in each direction is stored on the cell faces, whereas gas density and hydrodynamic pressure are stored at the cell centre. In order to compute the energy spectrum using (4.1), the gas momentum and velocity need to be computed in all directions at the cell centre. This requires a couple of interpolations. In the following, the process to interpolate the gas velocity field from the cell faces to cell centres is explained: (i) interpolate density from cell centres to cell faces using linear interpolation; (ii) compute velocity at cell faces,
$u_{j}=(\unicode[STIX]{x1D70C}u)_{j}/\unicode[STIX]{x1D70C}_{face_{j}}$
; (iii) take the three-dimensional Fourier transform of each velocity component,
$\hat{u} _{j}={\mathcal{F}}(u_{j})$
; (iv) shift
$\hat{u} _{j}$
by half a grid size,
$\hat{u} _{j}^{shifted}=\exp (-\text{i}\unicode[STIX]{x0394}x_{j}/2)\hat{u} _{j}$
; and (v) take the inverse Fourier transform of
$\hat{u} _{j}^{shifted}$
to obtain
$u_{j}$
at the cell centre.
Note that the momentum interpolation is similar, except that the first two steps are not required. The advantage of the above interpolation scheme is that the divergence-free condition in the incompressible limit is preserved. Therefore, in the limit, the Fourier transform of velocity is perpendicular to the modified wavenumber vector. Note that here a central difference is used to calculate the divergence. Since the Fourier transform of the velocity is unchanged during the interpolation (it is only shifted), the orthogonality condition will be preserved. For a general flow with variable density, the scheme described above does not produce artificial divergence.
Particle-to-gas heat transfer introduces turbulent fluctuations at high wavenumbers through flow dilatation. To study the turbulence modification caused by heat transfer from particles, we investigate the evolution of the divergence-free modes of the flow separately. This is accomplished by decomposing the gas momentum and velocity as follows:
$$\begin{eqnarray}u_{i}=u_{i}^{\unicode[STIX]{x1D735}}+u_{i}^{\prime },\quad \unicode[STIX]{x1D70C}u_{i}=(\unicode[STIX]{x1D70C}u)_{i}^{\unicode[STIX]{x1D735}}+\unicode[STIX]{x1D70C}u_{i}^{\prime }.\end{eqnarray}$$
The superscript
$\unicode[STIX]{x1D735}$
denotes the divergence-free part of a field, i.e. a Helmholtz–Leray decomposition (Bladel Reference Bladel1958). Accordingly, the TKE spectrum due to the divergence-free modes can be defined as follows:
$$\begin{eqnarray}E_{k}^{\unicode[STIX]{x1D735}}=\iint \text{Re}\{{\mathcal{F}}^{\ast }\{(\unicode[STIX]{x1D70C}u)_{j}^{\unicode[STIX]{x1D735}}\}{\mathcal{F}}\{u_{j}^{\unicode[STIX]{x1D735}}\}\}.\end{eqnarray}$$
In the remainder of the paper, the solid lines denoting energy spectra refer to
$E_{k}$
, and the dashed lines refer to the contribution of the divergence-free part,
$E_{k}^{\unicode[STIX]{x1D735}}$
.
4.1 Energy spectra with particle heating
For all of the results presented in this paper, a
$256^{3}$
grid is used. With this resolution, the Kolmogorov length scale is resolved, and the Eulerian grid size is much larger than the particle diameter. The latter condition is required for a valid point particle approximation. In figure 3(a) the normalized turbulent energy spectra after five large-eddy turnover times is compared between two simulations, one with a
$256^{3}$
grid and one with
$512^{3}$
. As illustrated, the wavenumber for which the energy accumulation begins is independent of the grid resolution. This is also evident from figure 3(b). Note that the amount of energy accumulated in the divergence part of the velocity field is approximately equal in both cases. However, since the
$256^{3}$
case consists of only half of the wavenumbers compared to the
$512^{3}$
case, it has higher energy per wavenumber values at the tail of the spectrum.
Figure 3. Comparison of the results after five eddy turnover times using a
$256^{3}$
grid and a
$512^{3}$
grid: (a) the normalized TKE spectra; (b) the divergence-free fraction of TKE.
Figure 4. The normalized TKE spectra for a decaying HIT at every large-eddy turnover time: (a) no particle heating; (b) nominal particle heating; (c) fraction of the TKE due to divergence-free modes for the case with nominal particle heating,
$\unicode[STIX]{x1D6FC}=1$
. The arrow shows the direction in time.
In the absence of a source, the HIT decays in time. In figure 4(a) the normalized TKE spectra for a decaying HIT is plotted for 10 large-eddy turnover times. Turbulence attenuation occurs at all wavenumbers. Figure 4(b) shows the evolution of the normalized TKE spectra in a decaying HIT in the presence of particle heating. Particle–turbulence interaction results in preferential concentration of particles (Eaton & Fessler Reference Eaton and Fessler1994). As we shall see in the analysis of energy transfer, heated particles transfer energy to the gas phase inhomogeneously. This inhomogeneous expansion strengthens the turbulence at all wavenumbers, but more significantly at higher wavenumbers. While the overall turbulence decays, particle heating results in extra energy at high wavenumbers compared to the standard HIT without heating. This excess energy manifests itself as plateaus in the energy spectra in the high-wavenumber dissipative range as shown in figure 4(b) (denoted by solid lines). For the case shown, the plateau of energy in wave space is the representation of divergence created due to heat released at the location of particles. Interestingly, the wavenumber span of the plateau increases as time progresses, reflecting an accumulation effect over time. A similar behaviour was observed in turbulent combustion (Kolla et al. Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014) where the energy spectrum displayed a plateau at high wavenumbers corresponding to dilatation generated at the flame scales.
In figure 4(b) dashed lines show
$E_{k}^{\unicode[STIX]{x1D735}}$
at different times. In the case with no particle heating (figure 4
a), the energy in the divergence-free modes is equal to the total energy and hence the solid and dashed lines are coincident. However, with particle heating the contribution of the divergence-free modes decreases at higher wavenumbers as evident in figure 4(b). In other words, the particle-to-gas heat transfer modifies the dilatational modes of the flow more significantly. In figure 4(c) the ratio of the divergence-free modes of the TKE to the total TKE is plotted as a function of wavenumber at different times. These Z-shaped curves shift towards the left along the
$x$
-axis as time progresses, indicating that, while energy is increasingly accumulating in the dilatational modes, the wavenumber range over which this occurs continues to broaden.
Figure 5. The effect of total particle heat absorption on a decaying particle-laden HIT. Four cases are shown: no heating; particle heating with nominal
$N_{p}$
; particle heating with
$\times 4$
particles; and particle heating with
$\times 16$
particles. (a) The normalized total TKE as a function of time. (b) The normalized TKE spectra after five large-eddy turnover times. (c) The fraction of the TKE due to divergence-free modes after five large-eddy turnover times.
To study the effect of total heat absorption by particles, we consider three different cases. One way to increase the total heat absorption is by increasing the number of particles. This will also proportionally increase
$\unicode[STIX]{x1D6FC}$
, the ratio of eddy time to the time scale associated with thermal energy increase. We have considered cases with four and 16 times more particles compared to the nominal number of particles. In figure 5(a) the temporal evolution of the total kinetic energy is shown for different levels of heating. Evidently, as the total heat absorption increases, the mean TKE decays more slowly.
Figure 5(b,c) indicates that the rate of heating can significantly affect the nature of TKE modification. Interestingly, for cases with higher total heating rate (which correspond to
$\unicode[STIX]{x1D6FC}=4$
and 16), turbulence modification shifts from dilatational modes to divergence-free modes. This is unintuitive since the primary impact of heating is to introduce dilatational modes in the flow. At the highest heating rate studied, the modification of TKE in dilatational modes can be ignored compared to that in divergence-free modes. This observation indicates that the energy transfer from divergence modes to divergence-free modes is highly nonlinear and becomes more and more efficient as the total heating magnitude is increased.
The generation of turbulence due to inhomogeneous heating results also in higher dissipation rate at high wavenumbers. In figure 6 the TKE dissipation rate is shown after five large-eddy turnover times for different levels of particle heating, similar to figure 5. In all cases, the maximum dissipation rate is well captured. Similar to the TKE spectra, for cases with more total heat absorption, we observe greater dissipation at smaller scales.
Figure 6. TKE dissipation after five large-eddy turnover times for different levels of total heat absorption similar to figure 5.
Figure 7. The normalized TKE spectrum after five eddy turnover times for three cases: no particle heating;
$N_{p}=5.12\times 10^{6}$
and
${\mathcal{R}}=0.34$
; and
$N_{p}=1.28\times 10^{6}$
and
${\mathcal{R}}=1.36$
.
An alternative way to increase
$\unicode[STIX]{x1D6FC}$
is to increase the particle heating intensity
${\mathcal{R}}$
in (3.6), while holding the number of particles fixed. In figure 7 two different scenarios with the same
$\unicode[STIX]{x1D6FC}$
are compared with the ‘no particle heating’ case. In the first case the number of particles is increased from nominal value by a factor of 4, while the heating intensity is maintained at the nominal level. In the second case, instead, the particle heating intensity is increased from nominal value by a factor of 4, while the number of particles is maintained at the nominal value. All other parameters are kept as the nominal values. The results show similar qualitative impact on the energy spectra.
4.2 Particle preferential concentration effect
As time progresses, three different phenomena occur simultaneously: turbulence decays, more heat is absorbed and Stokes number decreases. As time progresses, turbulence decays and therefore the Kolmogorov time scale increases. This results in a smaller Stokes number. The initial Stokes number is
$St_{\unicode[STIX]{x1D702}}=5.3$
. Therefore, as time progresses, the Stokes number approaches unity and stronger clusters of particles form (Eaton & Fessler Reference Eaton and Fessler1994). In figure 8 snapshots of particle distribution are shown after one, three and five large-eddy turnover times for the case with the nominal number of particles and heating intensity. Figure 8 demonstrates higher preferential concentration and particle temperature at later times.
Figure 8. Snapshots of particle distribution in a tiny slice of the domain after one, three and five large-eddy turnover times for the case with nominal particle heating. Particles are coloured by their non-dimensionalized temperature.
There are various ways to quantify the level of preferential concentration, which shows how the particle distribution differs from that of a random process: Fessler, Kulick & Eaton (Reference Fessler, Kulick and Eaton1994) used a box counting method; Reveillon & Demoulin (Reference Reveillon and Demoulin2007) used the energy spectrum of a particle density field; Monchaux, Bourgoin & Cartellier (Reference Monchaux, Bourgoin and Cartellier2010) used a Voronoi analysis; Reeks (Reference Reeks2014) used the full Lagrangian method (FLM), which provides a measure of the deformation of a volume of particles; and Sundaram & Collins (Reference Sundaram and Collins1997) and Yoshimoto & Goto (Reference Yoshimoto and Goto2007) used radial distribution functions (RDFs). Furthermore, in some of the previous works using the above methods, a characteristic length scale for particle clusters is introduced. For example, Aliseda et al. (Reference Aliseda, Cartellier, Hainaux and Lasheras2002) suggest a cluster size of order 10 Kolmogorov length scales. However, an investigation by Monchaux et al. (Reference Monchaux, Bourgoin and Cartellier2010) implies that clusters do not have a single characteristic size.
Figure 9. (a) Particle RDF for the case with nominal total heat absorption at every large-eddy turnover time. (b) Particle RDF for different levels of total heat absorption after five large-eddy turnover times. (c) Spatial autocorrelation of the gas temperature for different levels of total heat absorption after five large-eddy turnover times.
In this study we use RDFs (or pair correlation) to quantify the level of preferential concentration. RDF indicates the probability of finding pairs of particles with a certain radial distance compared to that of a random process. In the review by Monchaux, Bourgoin & Cartellier (Reference Monchaux, Bourgoin and Cartellier2012), pair correlation is recognized as the most widely used method to identify the preferential concentration. In figure 9(a), RDFs for particles with nominal
$N_{p}$
are depicted at different times. Consistent with figure 8, we observe evolution of stronger preferential concentration with time. This is because, as the turbulence decays, the Stokes number becomes closer to unity.
In figure 9(b) we show RDFs for different levels of heat absorption. For smaller levels of heating, the RDF is not influenced by heating, since the effect of heating on turbulence is weak, as seen in figure 4. We observe that higher heating tends to shift the system towards smaller clusters. The smaller magnitude of the clustering effect can be partially attributed to the diverging flows acting to spread out clusters in particle-crowded zones which naturally receive significant heating. This observation indicates that high heating not only modifies turbulence, but also influences the way that turbulence modulates particle clusters as seen in variation of RDF. Furthermore, we computed the radial autocorrelation of the gas temperature field, shown in figure 9(c). We observe a similar length scale trend as that in the particle RDFs. For the cases with a higher amount of total heat absorption, we observe slightly smaller clusters of hot gas consistent with the trends observed in figure 9(b).
Similar to the TKE spectrum defined in (4.1), we can define the gas temperature one-dimensional spectrum as
$$\begin{eqnarray}{\mathcal{T}}_{k}=\iint {\mathcal{F}}^{\ast }\{T_{g}\}{\mathcal{F}}\{T_{g}\}.\end{eqnarray}$$
In figure 10(a) the temporal evolution of
${\mathcal{T}}_{k}$
for the nominal case is shown. Initially, there is no variation in gas temperature, and thus
${\mathcal{T}}(t=0)=0$
. Over time, due to heating,
${\mathcal{T}}$
increases at all wavenumbers, but the increase is more significant at small wavenumbers. Note that the rate of heat injection to each particle is kept fixed in time. This constrains the thermal gradients around particles (i.e. the smallest scales) to be fixed such that the incoming heat is balanced with thermal diffusion. This explains the collapse of
${\mathcal{T}}_{k}$
for different times at large wavenumbers. In figure 10(b)
${\mathcal{T}}_{k}$
is shown for various cases after five large-eddy turnover times, indicating that with the higher total particle heat absorption we observe larger variations in
$T_{g}$
. From figure 10(b) we observe that increasing the particle heating intensity compared to increasing the number of particles (while
$N_{p}{\mathcal{R}}$
remains the same) results in larger
${\mathcal{T}}_{k}$
for both small and large scales.
Figure 10. Non-dimensional one-dimensional gas temperature spectra: (a) nominal total heat absorption after one, two, three, four and five large-eddy turnover times; (b) different levels of total particle heat absorption after five large-eddy turnover times.
4.3 Influence of dilatation
It is evident from figures 4 and 5 that dilatation plays a central role in the modification of the energy spectra. As noted in § 2, the role of dilatation in the balance equation for the energy spectrum is encapsulated in the pressure-dilatation term. In figure 11 we present one-dimensional spectra at selected instances for different cases: no particle heating, particle heating with nominal
$N_{p}$
and
${\mathcal{R}}$
, particle heating with
$N_{p}\times 4$
,
$N_{p}\times 16$
,
${\mathcal{R}}\times 4$
and
${\mathcal{R}}\times 16$
. Figure 11(a,d) illustrates the normalized one-dimensional spectra of the pressure-dilatation term as in (2.2) after one and five large-eddy turnover times. In the absence of particle heating, flow is divergence-free, and therefore the pressure-dilatation term vanishes. Figure 11(a,d) shows that, in the presence of particle heating, kinetic energy is injected to the flow through pressure-dilatation at a range of scales that are consistent with the cluster length scales observed in figure 9. We will verify this conclusion more rigorously in § 5 via investigating different cluster regimes by systematically changing the particle Stokes number. We here briefly note two trends consistent with this drawn conclusion. First, the peak in the pressure-dilatation term also shifts slightly to the right from time 1 to 5. This is due to the decay of turbulence, and therefore decrease of particle Stokes number, which results in smaller particle cluster length scales as seen in figures 8 and 9(a). Second, the peak wavenumber in the pressure-dilatation spectrum shifts slightly to smaller length scales with higher amount of total particle heating. This is consistent with figure 9(b), indicating that higher heating results in smaller clusters.
4.4 Spectral energy transfer
Since the mechanism of turbulent energy injection at wavenumbers corresponding to length scales of the particle clusters is evident through the dilatation influence, it is interesting to study the energy transfer across wavenumbers. The energy cascade can be modified by the pressure-dilatation effect. This will have consequences, particularly from an LES modelling point of view. With this motivation, we compute the energy transfer term and its evolution in time. Note that, in the turbulent kinetic energy equation, the energy transfer term as shown in (2.3) is on the left-hand side, while the pressure-dilatation term is on the right-hand side. Therefore, we consider the negative of the term introduced in (2.3) to be consistent with pressure-dilatation results.
Figure 11(b,e) shows the one-dimensional spectra of the energy transfer term for different cases at the selected instances. For all cases, the energy transfer term crosses over from negative to positive values beyond a threshold wavenumber, demonstrating energy transfer from large scales to small scales. Therefore, for the cases studied in this work, the presence of particle heating does not alter the overall direction of energy transfer (from low to high wavenumbers). However, the crossover wavenumber shifts to smaller scales with increasing net particle heat absorption, consistent with figure 5(c). This trend is even more pronounced at the later instant of five large-eddy turnover times, indicating that the cumulative effects of particle heating modify the energy transfer. Furthermore, the total amount of energy transfer from large scales to small scales grows as the net particle heat absorption increases. Note that the results of figure 11 shows that the spectral behaviour of the flow is similar for cases with an identical total heat absorption (
${\sim}N_{p}{\mathcal{R}}$
).
Figure 11(c,f) shows the superposition of the pressure-dilatation and energy transfer terms. The superposition plots collapse at scales larger than the particle cluster size. Therefore, the presence of particle heating has a slight influence at scales larger than the particle cluster size, as we also observed in figure 5(b).
5 Effect of particle Stokes number
In this section, we consider new cases with different Stokes number to study the effect of Stokes number on the results presented in previous sections. The nominal Stokes number at the beginning of turbulence decay is
$St_{\unicode[STIX]{x1D702}}\simeq 5$
. We consider two other cases with 10 times larger and 10 times smaller Stokes number. While we change the Stokes number, we keep all other parameters constant (including the thermal Stokes number). Therefore, changing Stokes number only affects the spatial distribution of particles. For all cases in this section we consider the maximum heating (
$16\times$
the base heating).
Figure 12. Particle distribution for cases with initial
$St=0.5$
,
$St=5$
and
$St=50$
, respectively, from left to right after five large-eddy turnover times. Only particles within a tiny slice are presented.
Figure 13. (a) Radial distribution function; and (b) gas temperature spatial autocorrelation for cases with initial
$St=0.5$
,
$St=5$
and
$St=50$
after five large-eddy turnover times.
As is visually verified in figure 12, after five large-eddy turnover times only the case with initial Stokes number
$\simeq 5$
demonstrates a significant particle preferential concentration (note that Stokes numbers decrease with time). This can be quantified by plotting the RDF similar to figure 9(a). The RDF plots for different cases are shown in figure 13(a). The case with initial Stokes number
$\simeq 5$
has the maximum deviation (
${\sim}100\,\%$
) compared to the other cases (18 % and 3 % for the cases with initial Stokes number 0.5 and 50, respectively). Therefore, the case with initial Stokes number
$\simeq 5$
has the highest level of preferential concentration. Moreover, figure 13(a) demonstrates that the length scale of clusters increases with particle Stokes numbers. Figure 13(b) shows the spatial autocorrelation of the gas temperature field after five large-eddy turnover times. Consistent with the length scales observed in figure 13(a), the length scales associated with the gas temperature field also increase with Stokes number.
Figure 14 shows the turbulent kinetic energy spectra for cases with different Stokes number after five large-eddy turnover times. As expected, the case with initial Stokes number 5 (i.e. highest level of preferential concentration) experiences the most significant TKE enhancement. In addition, all cases have approximately the same energy at the highest wavenumber. This observation can be explained by the fact that, given that these wavenumbers are dominated by dilatational modes, they represent local thermal expansions around individual particles. Since all cases involve the same number of particles with the same heating rates, similar TKE magnitudes are expected at the high wavenumbers.
Furthermore, from figure 14 we observe that the modification of the divergence-free modes (dashed lines) to the total kinetic energy is maximum for the case with the highest level of preferential concentration.
Figure 14. The normalized TKE spectra after five large-eddy turnover times for cases with initial
$St=0.5$
,
$St=5$
and
$St=50$
. Solid and dashed lines represent the total and divergence-free parts, respectively.
In figure 15(a,b) we present the pressure-dilatation and energy transfer terms for all three cases after five large-eddy turnover times. The range of length scales in which the pressure-dilatation term reaches its maximum increases with Stokes number. This is consistent with the trend of particle cluster size observed in figure 13. In addition, from figure 15(b) we observe that the deviation of the energy transfer term from that of the no-heating case is maximum for the case with the highest level of preferential concentration (i.e. the case with initial Stokes number 5). We also observe that the largest shift in the cross-over frequency occurs for the case with initial Stokes number 5.
Figure 15. (a) Normalized pressure-dilatation term; and (b) negative of energy transfer term for cases with initial
$St=0.5$
,
$St=5$
and
$St=50$
after five large-eddy turnover times.
6 Summary and discussion
Heating of the particles imparts dilatation to the carrier fluid, which introduces energy into the dilatational modes in a flow that is otherwise divergence-free. For systems with high heating rate, characterized by the parameter
$\unicode[STIX]{x1D6FC}$
(ratio of eddy time to thermal energy increase time), heating affects both the divergence and divergence-free modes of turbulence. Both the pressure-dilatation term, which encapsulates the dilatation source, and the energy transfer term increase in magnitude with increase in heating rate parameter
$\unicode[STIX]{x1D6FC}$
.
We found that the influence of particle heating, and ensuing heat transfer, on the energy cascade of particle-laden turbulent flows, interestingly, is analogous to that of heat release in turbulent flames. In this respect, the case with Stokes number of
$O(1)$
seems most relevant, since this scenario represents strong particle clustering with filament-like zones (see figure 8) of high particle concentration. Specifically, accumulation of energy at high wavenumbers was similarly observed in turbulent combustion where the energy spectrum displayed a plateau at wavenumbers corresponding to the flame scales (Kolla et al.
Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014). Hence, heated particle-laden turbulence may represent a surrogate model to study the spectral evolution of turbulence in reacting flows.
A surprising outcome of this investigation was revealed by decomposition of TKE spectra in terms of dilatational and divergence-free modes. We found that, for low heating rate, TKE modification is more pronounced in dilatational modes. Higher heating, to the level that the heating rate becomes faster than the eddy turnover time, results in significant modification of divergence-free modes. Under such heating rates, only a small fraction of TKE modification involves divergence modes.
The energy transfer analysis presented in § 4.4 shows significant modifications in the TKE spectra when particles are heated. As a result, backscatter of energy is possible regardless of whether decaying or sustained turbulence is concerned. The presented results can ultimately help with the development of SGS models for LES calculations, in which a wide range of flow structures affecting particles and thus the heat release zones are unresolved. In this work, we quantitatively presented the effects of heating in the most extreme clustering conditions. Future studies should focus on the quantification of these effects in a general
$St$
–
$\mathit{Re}$
map for relevant regimes of
$\unicode[STIX]{x1D70F}_{th}/\unicode[STIX]{x1D70F}_{p}$
. Based on such maps, one can identify the needs for model development in an LES setting. Additional work is needed to investigate how these models are applicable to scenarios with turbulence inhomogeneity such as the turbulent channel flow.
A possible extension of this work is to relax some of the simplifying assumptions. For instance, one could consider higher mass loading ratios and two-way momentum coupling as discussed in Elghobashi & Truesdell (Reference Elghobashi and Truesdell1993) and Horwitz & Mani (Reference Horwitz and Mani2016). Moreover, the effects of temperature-dependent gas properties (e.g. dynamic viscosity) on the conclusions presented in this study could be examined. Under these scenarios, one may need to reassess the contribution of all terms in the energy balance equation as presented in Kolla et al. (Reference Kolla, Hawkes, Kerstein, Swaminathan and Chen2014). Additionally, using the Helmholtz decomposition (see Bladel Reference Bladel1958) of the velocity field, the effects of divergence and divergence-free modes on the dynamics of energy transfer can be investigated separately.
In this study, we investigated how heating of preferentially concentrated particles brings about turbulence spectra modification. However, the reverse effect of such modification on the formation and strength of particle clustering (as can be inferred by figure 9 b) remains to be studied in future works.
Acknowledgements
We would like to acknowledge the anonymous referee for comments regarding the effects of particle Stokes number on TKE modification. This work was supported by the US Department of Energy under the Predictive Science Academic Alliance Program 2 (PSAAP2) at Stanford University. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the US Department of Energy under contract DE-AC04-94AL85000. The work at Sandia National Laboratories was supported by the Division of Chemical Sciences, Geosciences and Biosciences, Office of Basic Energy Sciences of the US Department of Energy.
