3.1. Characteristic time scales and non-dimensional numbers

The characteristic response time of the particle temperature to the surrounding fluid temperature field is given by the particle thermal relaxation time $\tau _{\theta }$ in (2.3). On the other hand, the characteristic response time of the fluid temperature field to the temperature of suspended particles is given by the fluid thermal relaxation time,

(3.1)\begin{equation} \tau_{\theta f} = (4 {\rm \pi}\kappa n_p r)^{{-}1},\end{equation}

where $n_p$ is the particle number density ($\tau _{\theta f}$ is referred to as the ‘gas thermal relaxation time’ in Pouransari & Mani (Reference Pouransari and Mani2018)). Here, $\tau _{\theta f}$ is inversely proportional to the particle number density $n_p$ and the particle radius $r$. This formula is similar to that for the phase relaxation time, an important time scale in cloud physics that characterizes the particle–turbulence interaction by condensation/evaporation (Politovich & Cooper Reference Politovich and Cooper1988; Korolev & Mazin Reference Korolev and Mazin2003; Kostinski Reference Kostinski2009). Recently, a time scale of similar form was found for the particle–turbulence interaction by the drag force (Saito et al. Reference Saito, Watanabe and Gotoh2019b).

There are three important non-dimensional numbers that characterize the two-way thermal interaction between particles and fluid. First, the particle-to-fluid heat capacity ratio is expressed as the ratio of $\tau _{\theta }$ to $\tau _{\theta f}$ as

(3.2)\begin{equation} \phi_{\theta} = \frac{C_p}{C_0} = \frac{c_p \phi_m}{c_0} = \frac{\tau_{\theta}}{\tau_{\theta f}}, \end{equation}

where $C_0$ and $C_p$ are the heat capacities of fluid and particles in the domain, respectively, $\phi _m$ is the mass-loading parameter given as

(3.3)\begin{equation} \phi_m = \frac{\rho_p \phi_v}{\rho_0} = \frac{\rho_p n_p (4 {\rm \pi}r^{3} /3)}{\rho_0}, \end{equation}

and $\phi _v$ is the volume fraction. Second, the ratio of $T$ to $\tau _{\theta f}$, where $T$ is the large-eddy turnover time for the turbulent velocity field, is represented by the Damköhler number as

(3.4)\begin{equation} Da = \frac{T}{\tau_{\theta f}}. \end{equation}

Here, $T={\mathcal {L}}/u_{rms}$ where ${\mathcal {L}}$ is the integral scale and $u_{rms}$ is the root-mean-square velocity (see Appendix D for definitions of other turbulence parameters). Third, the ratio of $\tau _{\theta }$ to $T$ is represented by the thermal Stokes number for the large-scale flow as

(3.5)\begin{equation} St_{L} = \frac{\tau_{\theta}}{T}, \end{equation}

where the subscript $L$ indicates that $St_{L}$ is that for the large-scale eddies of turbulent flow. From (3.4) and (3.5), $\phi _{\theta }$ is also expressed by $Da$ and $St_L$ as

(3.6)\begin{equation} \phi_{\theta} = Da St_L. \end{equation}

Note that the large-eddy turnover time, $T$, is used in the definitions of (3.4) and (3.5) because the largest eddies are considered to be the most important in the thermal coupling. This conjecture was first proposed for the particle–turbulence interaction by condensation/evaporation and was confirmed by DNSs of cloud turbulence (Lanotte, Seminara & Toschi Reference Lanotte, Seminara and Toschi2009; Sardina et al. Reference Sardina, Picano, Brandt and Caballero2015; Götzfried et al. Reference Götzfried, Kumar, Shaw and Schumacher2017; Kumar et al. Reference Kumar, Götzfried, Suresh, Schumacher and Shaw2018; Saito et al. Reference Saito, Gotoh and Watanabe2019a). Later this assumption was shown to also hold for the particle–turbulence interaction by the drag force (Saito et al. Reference Saito, Watanabe and Gotoh2019b). Note also that the importance of the Damköhler number was recently demonstrated by Pouransari & Mani (Reference Pouransari and Mani2018) for the particle-to-fluid heat transfer in particle-laden turbulence. They referred to the inverse of the Damköhler number ($Da^{-1} = \tau _{\theta f}/T$) as the heat-mixing parameter, as opposed to the commonly used heat-mixing parameter defined using the particle thermal relaxation time $\tau _{\theta }$. In the present study, we use the term ‘Damköhler number’ for (3.4) in order to emphasize the relation with cloud-turbulence research. Originally, the Damköhler number came from chemical engineering and described the ratio of the turbulent flow time scale and the chemical reaction time scale (Dimotakis Reference Dimotakis2005). In the cloud physics literature, a similar parameter was introduced to examine the role of cloud entrainment and turbulent mixing in the growth of cloud droplets by condensation/evaporation. Since the growth of cloud droplets has some similarities to a chemical reaction, the Damköhler number term has been used (Shaw Reference Shaw2003; Devenish et al. Reference Devenish2012). Saito et al. (Reference Saito, Watanabe and Gotoh2019b) introduced another similar parameter for momentum coupling and referred to it also as the Damköhler number. In order to distinguish a variety of Damköhler numbers, it is useful to add an additional adjective before it, for example: the supersaturation Damköhler number for cloud physics (as used by Shaw Reference Shaw2003); the inertial Damköhler number for momentum coupling; the thermal Damköhler number for (3.4). In the present study, however, we simply refer to (3.4) as the Damköhler number, for brevity.

3.2. Modulation of fluid temperature fluctuations by particles

We apply the Langevin model in order to derive the analytical expression for modulation of fluid temperature fluctuations by particles. We assume that turbulence is statistically homogeneous, is in a statistically steady state and is excited by the external forces $\boldsymbol {f}_{\boldsymbol {u}}$ and $f_{\theta }$ in (2.4) and (2.6), respectively, that are Gaussian random, white in time and homogeneous isotropic. We also assume that, both in cases with and without particles, the injection rate of the temperature-fluctuation variance is kept at $\chi _0$.

For cloud turbulence, the Langevin equation models the time variation of the supersaturation along the Lagrangian trajectory of cloud droplets. The feedback effect from cloud droplets on the supersaturation field is expressed by the term with the phase relaxation time $\tau _c$, and the turbulent mixing is expressed by the term with the large-eddy turnover time $T$. (For more details, readers are referred to, for example, Chandrakar et al. (Reference Chandrakar, Cantrell, Chang, Ciochetto, Niedermeier, Ovchinnikov, Shaw and Yang2016).) In the present system, i.e. (2.1) to (2.7), the fluid temperature $\theta$ corresponds to the supersaturation, and the fluid thermal relaxation time $\tau _{\theta f}$ corresponds to $\tau _c$. Therefore, the Langevin model can be straightforwardly applied to the present system.

The situation under consideration is as follows. The number density of particles is $n_p$ and particles are uniformly randomly distributed in the domain. As described in Saito et al. (Reference Saito, Watanabe and Gotoh2019b), we divide the domain into small subdomains and assume that each subdomain includes only a single particle (Tagawa et al. Reference Tagawa, Mercado, Prakash, Calzavarini, Sun and Lohse2012). The volume of the subdomain can be estimated approximately as $1/n_p$. The particle number density $n_p$ is assumed to be sufficiently large so that the characteristic length of the subdomain is smaller than the smallest length scale of the fluid temperature fluctuation (the Batchelor scale $\eta _{B}$) and that the fluid temperature in the subdomain can be regarded as uniform. In the subdomain, the fluid temperature is $\theta$ and the particle temperature is $\theta _p$. As time goes on, the fluid in the subdomain is advected by the turbulent velocity field as a small Lagrangian fluid parcel. Since we neglect the momentum inertia of the particle, the particle moves in the same way as the surrounding fluid parcel. From (2.6), the evolution equation for the temperature of the Lagrangian fluid parcel, $\theta$, is modelled by the Langevin equation as follows:

(3.7)\begin{equation} \theta(t+{\textrm{d}}t) - \theta(t) ={-} \alpha c \left( \frac{\theta - \theta_p}{\tau_{\theta f}} \right){\textrm{d}}t - \alpha \left( \frac{\theta}{T} \right){\textrm{d}}t + \alpha \left( \frac{2 \theta_0^{2} \,{\textrm{d}}t}{T} \right)^{1/2} \eta(t), \end{equation}

where $\theta _0$ is the standard deviation of the fluid temperature fluctuation without the effect of particles, $\eta (t)$ is a Gaussian random number with zero mean and unit variance which is statistically independent of $\theta$ and $\theta _p$, and $\alpha$ and $c$ are constants that are determined depending on the actual situation of the turbulence. The equation for the particle temperature is the same as (2.2), i.e.

(3.8)\begin{equation} \frac{{\textrm{d}} \theta_p}{{\textrm{d}}t} = \frac{\theta - \theta_p}{\tau_{\theta}}. \end{equation}

Since (3.7) and (3.8) are linear in $\theta$ and $\theta _p$, we can derive the analytical expressions for the statistically steady state.

We first consider the particle temperature $\theta _p$. Multiplying (3.8) by $\theta _p$ and taking an average, we have

(3.9)\begin{equation} \frac{{\textrm{d}}}{{\textrm{d}}t}\left( \frac{1}{2} \langle \theta_p^{2} \rangle_{p} \right) = \frac{1}{\tau_{\theta}} ( \langle \theta \theta_p \rangle_{p} - \langle \theta_p^{2} \rangle_{p}), \end{equation}

where $\langle \, \rangle _{p}$ indicates an average over particles. Since the left-hand side is zero in a statistically steady state, we have

(3.10)\begin{equation} \langle \theta \theta_p \rangle_{p} = \langle \theta_p^{2} \rangle_{p}. \end{equation}

Multiplying (3.7) by $\theta _p$ and taking an average $\langle \, \rangle _{p}$, we have

(3.11)\begin{equation} \left\langle \theta_p \frac{{\textrm{d}} \theta}{{\textrm{d}}t} \right\rangle_{p} ={-} \alpha \left( \frac{ \langle \theta \theta_p \rangle_{p} }{T} \right), \end{equation}

since the feedback terms from particles cancel out due to (3.10) and $\langle \eta (t) \theta _p \rangle _{p} =0$. Multiplying (3.8) by $\theta$ and taking an average $\langle \, \rangle _{p}$, we have

(3.12)\begin{equation} \left\langle \theta \frac{{\textrm{d}} \theta_p}{{\textrm{d}}t} \right\rangle_{p} = \frac{ \langle \theta^{2} \rangle - \langle \theta \theta_p \rangle_{p} }{\tau_{\theta}}. \end{equation}

Summing (3.11) and (3.12), we obtain

(3.13)\begin{equation} \frac{{\textrm{d}}}{{\textrm{d}}t}\left( \frac{1}{2} \langle \theta \theta_p \rangle_{p} \right) ={-} \alpha \left( \frac{ \langle \theta \theta_p \rangle_{p} }{T} \right) + \frac{1}{\tau_{\theta}} (\langle \theta^{2} \rangle - \langle \theta \theta_p \rangle_{p} ). \end{equation}

In a statistically steady state, we have

(3.14)\begin{equation} \langle \theta \theta_p \rangle_{p} = \left( \frac{\alpha}{T} + \frac{1}{\tau_{\theta}} \right)^{{-}1} \frac{\langle \theta^{2} \rangle_{p} }{\tau_{\theta}}. \end{equation}

Combining (3.10) and (3.14) and using the definition of (3.5), we finally obtain

(3.15)\begin{equation} \frac{ \langle \theta_p^{2} \rangle_{p} }{ \langle \theta^{2} \rangle_{p} } = \left( 1 + \alpha St_L \right)^{{-}1}. \end{equation}

This equation predicts the variance ratio of the fluctuations of particle and fluid temperatures.

We next consider the fluid temperature $\theta$. Here, we consider the dissipation of the temperature fluctuation variance. Since the injection rate is kept at a constant $\chi _0$ with or without particles, for the case without particles, we have

(3.16)\begin{equation} \chi_0 \sim \frac{\alpha}{T} \theta_0^{2} \end{equation}

and, for the case with particles, we have

(3.17)\begin{equation} \chi_0 \sim \frac{\alpha}{T} \langle \theta^{2} \rangle_{p} + \frac{\alpha c}{\tau_{\theta f}} ( \langle \theta^{2} \rangle_{p} - \langle \theta \theta_p \rangle_{p} ) \end{equation}

using (3.7). Equating (3.16) and (3.17) yields

(3.18)\begin{equation} \frac{\alpha}{T} \theta_0^{2} = \frac{\alpha}{T} \langle \theta^{2} \rangle_{p} + \frac{\alpha c}{\tau_{\theta f}} (\langle \theta^{2} \rangle_{p} - \langle \theta \theta_p \rangle_{p}). \end{equation}

Substituting (3.10) and (3.15) into the above equation and using the definitions of (3.2)–(3.6), we finally obtain

(3.19)\begin{equation} \frac{\langle \theta^{2} \rangle_{p}}{\theta_0^{2}} = \left( 1+ \frac{ c \alpha \phi_{\theta}}{1+ \alpha \phi_{\theta} Da^{{-}1}} \right)^{{-}1}. \end{equation}

This equation predicts modulation of the fluid temperature fluctuation variance by particles as a function of the Damköhler number.

We consider two limits of (3.19). First, for the limit $Da \rightarrow 0$ with fixed $\phi _{\theta }$, we have

(3.20)\begin{equation} \frac{\langle \theta^{2} \rangle_{p}}{\theta_0^{2}} \rightarrow \left( 1 + c Da \right)^{{-}1} \quad (Da \rightarrow 0). \end{equation}

Second, for the limit $Da \rightarrow \infty$ with fixed $\phi _{\theta }$, we have

(3.21)\begin{equation} \frac{\langle \theta^{2} \rangle_{p}}{\theta_0^{2}} \rightarrow \left( 1 + c \alpha \phi_{\theta} \right)^{{-}1} \quad (Da \rightarrow \infty), \end{equation}

which indicates that modulation of the variance approaches a constant depending on $\phi _{\theta }$.

In summary, in a statistically steady state, (3.15) predicts the variance ratio of the particle and fluid temperature fluctuations, and (3.19) predicts modulation of the fluid temperature fluctuation variance by particles.

3.2.1. Modulation by particles with fixed temperature

An important limit case of the Langevin model given by (3.7) and (3.8) is the case in which the particle temperature fluctuation $\theta _p$ is fixed at zero. This physically corresponds to particles with large thermal inertia ($St_L \rightarrow \infty$). In this case, the Langevin model reduces to

(3.22)\begin{equation} \theta(t+{\textrm{d}}t) - \theta(t) ={-} \alpha c \left( \frac{\theta}{\tau_{\theta f}} \right){\textrm{d}}t - \alpha \left( \frac{\theta}{T} \right){\textrm{d}}t + \alpha \left( \frac{2 \theta_0^{2} \,{\textrm{d}}t}{T} \right)^{1/2} \eta(t),\end{equation}

which is equivalent to the Langevin model used in cloud turbulence research (except for constants $c$ and $\alpha$ – see (2) in Chandrakar et al. (Reference Chandrakar, Cantrell, Chang, Ciochetto, Niedermeier, Ovchinnikov, Shaw and Yang2016)). In a statistically steady state, we obtain the same result as (3.20), which is the same formula as (3) in Chandrakar et al. (Reference Chandrakar, Cantrell, Chang, Ciochetto, Niedermeier, Ovchinnikov, Shaw and Yang2016).

3.2.2. Particle temperature fluctuation for the one-way coupling case

Another important limit case of the Langevin model (3.7) and (3.8) is the case without the feedback effect of particles on the fluid temperature field, namely, the one-way coupling case. In this case, (3.7) reduces to

(3.23)\begin{equation} \theta(t+{\textrm{d}}t) - \theta(t) ={-} \alpha \left( \frac{\theta}{T} \right){\textrm{d}}t + \alpha \left( \frac{2 \theta_0^{2} \,{\textrm{d}}t}{T} \right)^{1/2} \eta(t), \end{equation}

and the equation for the particle temperature is the same as (3.8), as follows:

(3.24)\begin{equation} \frac{{\textrm{d}} \theta_p}{{\textrm{d}}t} = \frac{\theta - \theta_p}{\tau_{\theta}}.\end{equation}

In a statistically steady state, we obtain the same result as (3.15) for the particle temperature, as follows:

(3.25)\begin{equation} \frac{ \langle \theta_p^{2} \rangle_{p} }{ \langle \theta^{2} \rangle_{p} } = \left( 1 + \alpha St_L \right)^{{-}1}, \end{equation}

whereas $\langle \theta ^{2} \rangle _{p} = \theta _0^{2}$ for the fluid temperature because the coupling is one-way. Note that a relationship similar to (3.25) was recently derived by Carbone et al. (Reference Carbone, Bragg and Iovieno2019) (see their (4.5)).

4.1. Simulation set-up

In order to confirm the theoretical prediction obtained in the previous section, we conduct DNS of particles in turbulence using the cloud microphysics simulator, which is a DNS model developed for the large-scale computation of particle-laden turbulent flow (Gotoh, Suehiro & Saito Reference Gotoh, Suehiro and Saito2016; Saito & Gotoh Reference Saito and Gotoh2018). The numerical domain is a periodic cubic box with a length $L_{box}$ per side. Particles are regarded as spherical point particles. The fluid velocity and temperature at the particle position, $\boldsymbol {u} (\boldsymbol {x}_{pj} )$ and $\theta (\boldsymbol {x}_{pj})$, are calculated from the velocity and temperature fields at the surrounding eight grid points by linear interpolation, respectively. The same weighting is used to project the particle temperature, $\theta _{pj}$, onto the surrounding eight grid points to calculate the term $F_p$ in (2.7), i.e. the particle-in-cell method (Sundaram & Collins Reference Sundaram and Collins1996).

We numerically integrate the evolution equations (2.1), (2.2), (2.4) and (2.6) using the second-order Runge–Kutta scheme and the pseudospectral method for temporal and spatial discretization, respectively. The external forces $\boldsymbol {f}_{\boldsymbol {u}}$ and $f_{\theta }$ in (2.4) and (2.6) are homogeneous isotropic, applied to all wavenumber components within a shell of $1 \leq k L_{box} \leq 4$, and are produced by Gaussian random, white-in-time processes (Gotoh, Fukuyama & Nakano Reference Gotoh, Fukuyama and Nakano2002). The details of the external forces are provided in Appendix A. These simple forces are enough for the fundamental study. After a statistically steady state is achieved for each simulation, turbulence and particle statistics are calculated.

Numerical and physical parameters in the present DNS are summarized in table 1. We studied five values of the total grid number ($N_{grid}^{3}=32^{3}$–$512^{3}$) and 10 values of the particle radius ($r=1.0$–$1000\,\mathrm {\mu }$m). The total number of particles $N_p$ is set equal to the total grid number $N_{grid}^{3}$, corresponding to the particle number density $n_p = 125$ cm$^{-3}$. For the particle-to-fluid heat capacity ratio, we chose three values ($\phi _{\theta }=1.0$, $0.3$ and $0.1$) and conducted two additional simulations. In one simulation, the particle temperature was fixed at zero ($\theta _{pj}=0$) and, in the other simulation, the thermal coupling was one-way (the feedback from particles onto the fluid was neglected). In total, we conducted $5 \times 10 \times (3+2)=250$ kinds of simulations. Table 2 summarizes the particle radius $r$, the fluid thermal relaxation time $\tau _{\theta f}$, the volume fraction $\phi _v$ and the particle-to-fluid density ratio $\rho _p/\rho _0$ for each simulation.

Table 1.

Various numerical and physical parameters. The time step is $\Delta t = 1.0 \times 10^{-3}$ (s) for simulations with $N_{grid}=32$, $64$, $128$, and $256$ and is $\Delta t = 8.0 \times 10^{-4}$ (s) for simulations with $N_{grid}=512$. The box length is $L_{box}=N_{grid} \Delta x$. The total number of particles is $N_p= n_p L_{box}^{3} =N_{grid}^{3}$. The density of particles is $\rho _p=(3 \rho _0 c_0 \phi _{\theta })/(4 {\rm \pi}n_p c_p r^{3})$. For parameter $\phi _{\theta }$, ‘1.0 (one-way)’ indicates simulations with $\phi _{\theta }=1$ and one-way thermal coupling, and ‘$\theta _p$ zero’ indicates simulations with the particle temperature fixed at zero.

Table 2.

Particle parameters. Here, $r$ is the particle radius, $\tau _{\theta f}$ is the fluid thermal relaxation time and $\phi _v$ is the volume fraction. The particle-to-fluid density ratio ($\rho _p/\rho _0$, as determined by (3.2) and (3.3)) is listed for simulations with $\phi _{\theta }=1.0$, $0.3$ and $0.1$. The particle thermal relaxation time is given by $\tau _{\theta }=\phi _{\theta } \tau _{\theta f}$ ($\phi _{\theta }=1.0, 0.3, 0.1$).

In table 1, the values of $c_0$ and $c_p$ are taken from Saito et al. (Reference Saito, Gotoh and Watanabe2019a) and Sato, Deutsch & Simonin (Reference Sato, Deutsch and Simonin1998), respectively. From (3.2), the mass-loading parameters corresponding to $\phi _{\theta }=1.0$, $0.3$ and $0.1$ are $\phi _m = 2.6$, $0.78$ and $0.26$, respectively. Although $\phi _m$ is an important parameter for turbulence modulation by particles due to the drag force (Squires & Eaton Reference Squires and Eaton1993; Saito et al. Reference Saito, Watanabe and Gotoh2019b), $\phi _m$ is not important in the present study because momentum coupling is neglected and the particles move as fluid particles.

The important premise in the Langevin model in the previous section is that the largest eddies play the central role in the particle–turbulence interaction by the two-way thermal coupling. We check the validity of this premise using the same strategy as that used by Lanotte et al. (Reference Lanotte, Seminara and Toschi2009) (see also Saito et al. (Reference Saito, Watanabe and Gotoh2019b)). Namely, we fix the turbulence statistics for the smallest eddies and expand the number of grids, the domain size and the sizes of the largest eddies in order to investigate the interaction between the largest eddies and particles. Turbulence statistics without particles are summarized in table 3. Note that fluid temperature statistics are modulated by the effect of particles.

Table 3.

Mean turbulence parameters for cases without particles. Here, $N_{grid}$ is the number of grid points in physical space for each direction, $R_{\lambda }$ is the Taylor microscale Reynolds number, $u_{rms}$ is the root-mean-square velocity, $\epsilon$ is the mean energy dissipation rate per unit mass, $\mathcal {L}$ is the integral scale, $\lambda$ is the Taylor microscale, $\eta$ is the Kolmogorov length, $k_{max}\eta$ is the cutoff wavenumber normalized by the Kolmogorov length, $T$ is the large-eddy turnover time, $\tau$ is the Kolmogorov time, $\theta _{rms}$ is the root-mean-square temperature fluctuation ($=\theta _0$ in (3.7)) and $\chi$ is the mean dissipation rate for the temperature variance. See Appendix D for definitions of turbulence parameters.

In the present study, the system is characterized by five independent non-dimensional numbers: the Taylor microscale Reynolds number $Re_{\lambda }$; the Prandtl number $Pr=\nu /\kappa$ which is unity; the particle-to-fluid heat capacity ratio $\phi _{\theta } = \tau _{\theta }/\tau _{\theta f}$; the Damköhler number $Da = T/\tau _{\theta f}$ (or the thermal Stokes number for lerge scale flow $St_L = \tau _{\theta }/T$); and the thermal Stokes number defined by the Kolmogorov time, $St_{\theta }=\tau _{\theta }/\tau$, where $\tau$ is the Kolmogorov time. Either $Da$ or $St_L$ can be removed by (3.6).

For all of the following simulation results, we confirmed that $\langle \theta ^{2} \rangle _p$ is almost identical to $\langle \theta ^{2} \rangle _V$, where $\langle \, \rangle _V$ indicates the volume average (agreement is better for simulations with larger $N_{grid}$ and $N_p$, ranging from 97.1 % for $N_{grid}=32$ to 99.6 % for $N_{grid}=512$). In other words, the variance of the fluid temperature fluctuation does not change whether the average is taken over particles or over grid points. This is because the total number of particles $N_p$ is sufficiently large (same as the total grid number) and particles are almost uniformly randomly distributed in the domain in the present simulation. We also confirmed that the probability density functions of both the fluid and particle temperature fluctuations are all Gaussian.