Statistical properties of thermally expandable particles in soft-turbulence Rayleigh-Bénard convection

The dynamics of inertial particles in Rayleigh-Bénard convection, where both particles and fluid exhibit thermal expansion, is studied using direct numerical simulations (DNS) in the soft-turbulence regime. We consider the effect of particles with a thermal expansion coefficient larger than that of the fluid, causing particles to become lighter than the fluid near the hot bottom plate and heavier than the fluid near the cold top plate. Because of the opposite directions of the net Archimedes’ force on particles and fluid, particles deposited at the plate now experience a relative force towards the bulk. The characteristic time for this motion towards the bulk to happen, quantified as the time particles spend inside the thermal boundary layers (BLs) at the plates, is shown to depend on the thermal response time, τT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \tau_{T}$\end{document}, and the thermal expansion coefficient of particles relative to that of the fluid, K=αp/αf\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ K = \alpha_{p}/\alpha_{f}$\end{document}. In particular, the residence time is constant for small thermal response times, τT≲1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \tau_{T} \lesssim 1$\end{document}, and increasing with τT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \tau_{T}$\end{document} for larger thermal response times, τT≳1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \tau_{T} \gtrsim 1$\end{document}. Also, the thermal BL residence time is increasing with decreasing K. A one-dimensional (1D) model is developed, where particles experience thermal inertia and their motion is purely dependent on the buoyancy force. Although the values do not match one-to-one, this highly simplified 1D model does predict a regime of a constant thermal BL residence time for smaller thermal response times and a regime of increasing residence time with τT\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$ \tau_{T}$\end{document} for larger response times, thus explaining the trends in the DNS data well.


Introduction
Inertial particles in thermally driven flows are abundant in both nature and technological applications. In nature typical examples are aerosols in the atmospheric boundary layer [1], the dynamics of droplets in clouds [2,3] or plankton in oceanic flows [4,5], while in technological applications one can think of spray combustion [6,7] or solar collectors [8]. Particles in flows occur in a wide range of densities; while plankton and algae in the ocean have a density close to that of the carrier fluid, droplets in clouds are in general much heavier than the surrounding fluid. When the particle density is different from the fluid density, inertia will cause particle trajectories to deviate from the fluid streamlines, resulting in a non-homogeneous distribution of particles in the flow [9][10][11][12]. When heat transfer Contribution to the Topical Issue "Flowing Matter, Problems and Applications", edited by Federico Toschi, Ignacio Pagonabarraga Mora, Nuno Araujo, Marcello Sega. a e-mail: k.m.j.alards@tue.nl b e-mail: F.Toschi@tue.nl between particles and fluid is not instantaneous, also thermal inertia plays a role. Thermal inertia takes into account the time particles need to adjust their internal temperature to that of the surrounding fluid, which is typically referred to as the thermal response time.
The effect of thermal inertia will be visible in the temperature statistics of the particles. The larger the thermal response time of particles, the more the temperature of particles will deviate from the underlying fluid temperature (at the particle position). When considering a dilute suspension (where particles are not expected to influence the fluid flow or temperature) where the size of particles is independent of temperature, thermal inertia will not influence the motion of the inertial particles. However, when the volume of the particles does depend on the temperature, thermal inertia can drastically change the trajectories of the inertial particles. For example, bubbles in boiling convection will grow in the warmer spots of the flow and shrink in the cooler spots [13][14][15], affecting their buoyancy and therefore changing the upward and downward motion of these bubbles. This behavior is not restricted to bubbles, but also, e.g., trajectories of (phase changing) oil droplets [16] or gel-like particles in non-isothermal flows are expected to be influenced by thermal inertia.
Here we conduct a numerical study on the dispersion of thermally and mechanically inertial particles in Rayleigh-Bénard convection (RBC), a fluid layer heated from below and cooled from above. The typical flow structure in such a RBC set-up is a large-scale circulation (LSC) of rising hot fluid and descending cold fluid [17]. We focus on the higher end of the "soft-turbulence" Rayleigh-Bénard regime, according to the classification of [18].
The temperature dependence of the particle size is included as thermal expansion. This method is not restricted to bubbles but can also deal with fluid-fluid systems or gellike particles in non-isothermal flows. Not only the thermal expansion of particles, but also that of the fluid is taken into account, such that the volume of both particles and fluid increases (linearly) with increasing temperature. In particular, we study particles with an average density equal to that of the fluid and a thermal expansion coefficient larger than that of the fluid. In this setting, particles become lighter than the fluid near the hot bottom plate and heavier than the fluid near the cold top plate. This is expected to induce an enhanced upward or downward motion to the particles, respectively, on top of the motion of plumes near the plates. These plumes were shown to be able to transport inertial particles away from the plates, however (without thermal expansion) particles were eventually deposited at the plates again due to the gravitational force [19]. By including thermal expansion we expect particles to be transported towards the plates by the LSC, be deposited on the plate due to their mechanical inertia, and stay there for some characteristic residence time to then re-suspend due to their enhanced thermal expansion compared to the fluid.
In this way, thermal expansion of particles prevents them from definitively settling at the horizontal plates. In experiments, even a very small mismatch between fluid and particle density leads to particles getting deposited at the top and bottom plates (as, e.g., in [20]). This settling of particles will not only reduce the number of particles inside the bulk flow, but could also have significant effects on the heat transfer. The effect of thermally conductive particles with density very close to that of the fluid on the heat transfer in RBC was investigated experimentally by Joshi et al. [20]. Particles were found to settle at the walls, depleting the RBC bulk flow of particles and forming a porous layer at the plates that eventually would cause a decrease of the heat transfer. In numerical studies particles are often prevented from getting stuck at the plates by neglecting gravity [21][22][23], by pointing gravity in the direction parallel to the walls [24,25] or by removing particles from the flow as soon as they reach one of the plates [14,15,26]. Here, the larger thermal expansion coefficient of particles alone ensures that particles eventually move away from the plates again.
The dynamics of thermally inertial particles (without thermal expansion) in RBC has already been studied numerically in the limit of bubbles (light particles) [13][14][15] and in the limit of particles which are heavier than the fluid [26]. In these studies a two-way coupling approach is used, i.e. the feedback reaction of particles on the fluid velocity and temperature is included in the momentum and energy equations. It was found that these two-way coupled inertial particles significantly affect the heat transfer due to the mismatch between fluid and particle density. However, as a result of this density mismatch particles will get stuck at the horizontal plates. Here we consider particles with a temperature-dependent density but with an average density equal to that of the fluid. In this regime of density ratios particles are not expected to significantly influence the heat transfer and flow structures. A one-way coupling treatment is then sufficient. An example of a system where particles have a density very close to that of the fluid, but also a larger thermal expansion coefficient than the fluid, is a configuration of gel-like particles in water; these particles can consist of a rubber coating filled with a mineral or silicon gel [27]. We study how thermal inertia affects the dynamics and the distribution of such particles in RBC.
In the remainder of this paper we first introduce the numerical set-up in sect. 2.1 by explaining both the RBC flow set-up and the modeling of thermally and mechanically inertial particles. In sect. 3, we discuss our results, focusing on the distribution and dynamics of these thermally responsive particles in RBC. Results will be presented for a wide range of thermal response times and for different ratios between the thermal expansion coefficient of the particles and of the fluid. In the last sect. 4 we will summarize and conclude our findings.

Numerical methods
We study RBC, seeded with thermally and mechanically inertial particles, using direct numerical simulations (DNS). Below we will discuss both the numerical model for RBC and the modeling of the inertial particles.

Rayleigh-Bénard convection
In RBC a fluid is heated from below and cooled from above, inducing a buoyancy driven flow. Control parameters for the RBC set-up are the Rayleigh number, Ra = α f gΔT H 3 /(κν), and the Prandtl number, P r = ν/κ, with α f the thermal expansion coefficient of the fluid, g the gravitational acceleration, ΔT the temperature difference between the plates, H the height of the RBC cell and κ and ν the thermal diffusivity and the kinematic viscosity of the fluid, respectively. The numerical Rayleigh-Bénard set-up studied here is bounded by top and bottom walls and has periodic boundary conditions in the horizontal directions. The governing dimensionless equations are the incompressible Navier-Stokes and energy equations in the Boussinesq approximation: (1) grid points and to also ensure at least ten grid points in the thermal and viscous boundary layers (BLs), grid refinement is used in the vertical direction. The discretization is performed on a staggered grid using a second-order finite-difference scheme and for the integration a thirdorder Runge-Kutta method is applied. Details of the numerical scheme can be found in [28,29]. In this study the Prandtl number was chosen equal to P r = 6.7 (corresponding to water) and the Rayleigh number was chosen equal to Ra = 2 · 10 7 , i.e. at the higher end of the "softturbulence" regime, according to the classification introduce in [18]. All fluid properties, corresponding to the average fluid temperature T m = 0.5, are reported in detail in table 1.

Thermally expanding inertial particles
Particles which experience both thermal and mechanical inertia are evolved in the RBC flow. We treat these particles as point particles, a reasonable assumption when the radius of particles, r p , is smaller than the smallest length scale of the flow, η, the Kolmogorov length scale. Note that in RBC a second length scale is involved related to the temperature field; the Batchelor length η B . In the setup studied here this length scale is smaller than η, since η B = η/ √ P r ≈ 0.4η. To derive the equation for the thermal inertia it is additionally assumed that the thermal conductivity of the particles is much larger than that of the fluid, such that the Biot number of the particles is small, Bi 1, and temperature gradients within the particles can be neglected [30]. The equation for the velocity of one particle, u p , is based on the Maxey-Riley equation [31] and for the temperature of that particle, T p , the approach proposed by Michaelides in [30] is used, such that where u f (x p ) and T f (x p ) are the fluid velocity and the fluid temperature at the position of the particle, x p , respectively. Here β = ρ p /ρ f is the ratio between the density of that particle and the fluid density, Re p = 2r p |u p − u f (x p )|/ν is the particle Reynolds number and τ p and τ T are the viscous and thermal response times, respectively. These are defined as where γ = c p /c f is the ratio between the specific heats of the particle material, c p , and the fluid, c f , [30,31]. The forces, included on the right-hand side (rhs) of eq. (4), are the Stokes drag, the added mass (also responsible for the pre-factor on the left-hand side) and the gravitational force. In eq. (5), the term on the rhs is analogue to the drag force. Since the particles simulated here have a particle Reynolds number of about Re p ∼ 10 it is necessary to include non-linear effects in the drag forces, represented by the factors (1 + 0.15Re 0.687 p ) in eq. (4) [32] and (1 + 0.3Re [33]. The pressure gradient force and the Basset history force are not included in eq. (4), while these forces might be important in a system where particle and fluid density are similar and β ≈ 1 [34,35]. We verified that ignoring these terms is not influencing the (statistical) measures discussed in this paper and that the most important contributions actually come from the Stokes drag force, added mass force and the buoyancy force. For clarity we therefore choose to not include the Basset history force and the pressure gradient force. In the equation for the thermal inertia we neglect both the history force and the force analogue to the added mass contribution [30], again after verifying that the contribution of these terms is minor and that the most important contribution comes from the term analogue to the drag force.
As mentioned above eq. (5) is valid for Bi 1. The Biot number of the particles is defined as Bi = 2h p r p /k p , with h p and k p the heat transfer coefficient and the thermal conductivity of particles, respectively. It is possible to express h p in terms of a particle-Nusselt number, Nu p = 2r p h p /k f , with k f the thermal conductivity of the fluid. The Biot number thus becomes Bi = Nu p k f /k p , proportional to the ratio between the thermal conductivity of the fluid and the particles. In a suspension of solid particles in a fluid, k p k f and the assumption of Bi 1 is indeed valid. However, in fluid-fluid systems Bi ∼ O(1) making temperature differences between the core and the surface of particles possible. We expect that this will not significantly effect our results and will at most result in an additional delay in the heat transfer between the particle and the surrounding fluid. This will lead to a larger "effective" thermal response time and since results are presented in a wide range of thermal response times we expect that our results are also applicable to the case of Bi ∼ O (1).
Particles and fluid both exhibit thermal expansion with different thermal expansion coefficients, where the thermal expansion coefficient of the particles is chosen to be larger than that of the fluid, such that α p > α f . The densities of the fluid and the particles are assumed to decrease linearly with the temperature fluctuations of the fluid (T f = T f − T m ) and the fluctuations in the particle temperature (T p = T p − T m ): where the densities of the particles and the fluid at the average temperature are set to unity, ρ p (T p = T m ) = ρ f (T f = T m ) = 1, without loss of generality. Now also the density ratio is temperature dependent, as Due to the thermal expansion also the size of the particles depends on the temperature fluctuations. Under the assumption that temperature fluctuations are small (as also assumed by the Boussinesq approximation for eq. (2)) and by using the Taylor expansion, the radius of particles follows as where r p is the temperature-dependent radius, while r p is the radius of particles at T p = T m and where higher order terms have been ignored.
Since the viscous and thermal response times depend on both the density ratio and the particle radius, they have to be updated accordingly such that Table 2. Particle properties of the thermally responsive particles (TRP), simulated in Rayleigh-Bénard convection. Three different simulations are performed with tracers (family 0) and thermally responsive particles (families 1-9), for three different ratios between the thermal expansion coefficient of the particles and that of the fluid, K = α p/αf . Here rp, τp and β are the particle radius, the drag response time and the ratio between the particle and fluid density at the mean temperature T m, respectively. The properties of the different particle families at the average particle and fluid temperature, are reported at the bottom of the table, where γ = c p/cf , with cp and c f the specific heat of the particles and the fluid, respectively, and τ T is the thermal response time. where τ p and τ T are the particle and thermal response times at T p = T f = T m , respectively, and we again neglect higher order terms. To complete the implementation of thermal expansion, the parameters β, τ p and τ T , in eqs. (4) and (5) have to be replaced by the temperature-dependent variables β, τ p and τ T , respectively. Also the particle Reynolds number is now based on the temperature-dependent radius r p . The typical time these thermally responsive particles spend at the plate in order to adjust their density enough to escape the BLs, is expected to depend on the ratio between the specific heats of the particle material and the fluid, γ. Therefore we study particles in a wide range of thermal response times, τ T . On top of this, we introduce a key parameter for this study, K = α p /α f , being the ratio between the thermal expansion coefficient of the particle and that of the fluid. Three different values of this parameter K are studied: K = 1.1, K = 2 and K = 10, as also reported in table 2. The applications mentioned in the introduction, gel-like like particles in water and oil-water configurations, would fall in the range of 1.1 K 2. Here, K = 10 is added to also study a more extreme case. For each value of K, ten different particle families are included in the simulation; one family consisting of passive tracers and nine families of thermally responsive particles, with 0.05 ≤ τ T ≤ 10 as reported in table 2. These thermal response times correspond to a range of 0.13 ≤ γ ≤ 26. In general γ ∼ O(1) for solid-fluid or fluid-fluid systems. To give an estimate of the corresponding thermal response times; in a range of 0.3 ≤ γ ≤ 3 the thermal response times would be 0.12 ≤ τ T ≤ 1.2. Here we again add extreme values of both smaller and larger γ to understand how the systems converges in the limit of very small and very large thermal response times. In this parameter range the density ratio varies between 0.96 < β < 1.04. Within this range of density ratios particles are not expected to influence the flow structures and the heat transfer and therefore a one-way coupling approach is sufficient. In total nine different particle families are simulated for 300 dimensionless time units, where the number of particles in each family is 1.6 · 10 5 . A detailed overview of the particle properties is given in table 2.

Spatial distribution of thermally expandable particles
We investigate the dynamics of thermally responsive inertial particles in Rayleigh-Bénard convection, where we include thermal expansion of both particles and of fluid. In particular, the thermal expansion coefficient of particles is larger than that of the fluid, so that particles react to the temperature fluctuations stronger. Since in RBC the temperature gradients are largest in the thermal BLs while the temperature in the bulk fluctuates around the average temperature [17,36], particles are expected to distribute differently in the bulk than in the thermal BLs when thermal expansion is included. To study the vertical distribution of particles, we compute the particle number density, n i , as a function of z. First the RBC cell is subdivided into 250 horizontal slabs of size Δz = 0.004H, with central vertical position z i . The number density in each slab is computed as the time averaged number of particles in the slab divided by the slab volume; Finally, this number density is normalized by the total number density N tot /V tot , where N tot = 1.6 · 10 5 (for each particle family) and V tot = HL x L y . In summary this means n i = Ni Vi / Ntot Vtot . In fig. 1, we show n i for the three different values of K: K = 1.1, K = 2 and K = 10 and different values of τ T between τ T = 0.05 and τ T = 10. As a reference the distribution of fluid tracers is also shown with gray lines with crosses. As expected, fluid tracers are distributed uniformly such that n i = 1. Note that these fluid tracers have no thermal and mechanical inertia (τ T = τ p = 0) and that they are therefore not affected by thermal expansion. The thermal BL thickness, δ T = 0.022H, is computed as the position of the maximum root-mean-square temperature and is indicated in fig. 1 by the vertical black lines. First, we observe that the number of particles inside the thermal BL is increasing with increasing thermal response times compared to the uniform distribution n i = 1. Particles with a larger thermal response time need more time to heat up (cool down) at the bottom (top) plate, hence there will be more particles close to the plates on average. Furthermore, when comparing the three different panels, it is observed that this number of particles at the plate is larger for lower values of K. Particles with a larger thermal expansion coefficient compared to that of the fluid react very strongly to temperature fluctuations and even a small temperature change can lead to a huge change in their mass density. Consequently, particles move away from the plates faster and the number of particles at the plates decreases. For K = 2 and K = 10 we observe a regime where n i < 1 for τ T 2 and τ T 4, respectively. Here particles escape the BLs so fast that there is a depletion of particles in the thermal BLs, compared to the average distribution n i = 1. A depletion in the BLs results in an increase of particles in the bulk, indicated by the peaks in figs. 1(b) and (c) for z δ T , which become more prominent for larger K and for smaller τ T .
The particle number density, as shown in fig. 1, is an average quantity and does not give information on the particle distribution in the horizontal directions. To understand how particles distribute horizontally with respect to the typical temperature profiles in the RBC cell, we visualize the temperature field at z = 0.012H without particles in fig. 2(a) and with different types of thermally responsive particles with vertical position z p < 0.015H in figs. 2(b)-2(g), where particles are colored by their temperature. For each value of K (K = 1.1, K = 2 and K = 10), a situation with a low thermal response time of τ T = 0.1 and a situation with a large thermal response time of τ T = 4 are shown. First, when focusing on the effect of the thermal response time, it is observed that there are more particles at the plate for larger thermal response times (as already discussed above) and that particles with a lower thermal response time are only found in the colder spots. These particles have a temperature very close to that of the fluid and it is expected that colder heavier particles stay at the plates longer, explaining why in this regime colder particles are found, clustered in the colder spots of the fluid at the bottom plate in panels (b), (d) and (f). For larger thermal response times, particle and fluid temperature are less correlated and, especially for lower values of K, particles are less restricted to the colder areas of the fluid, see panels (c), (e) and (g).

Temperature statistics
From fig. 2 we expect that the distribution of particles is related to the temperature of the particles, relative to the temperature of the surrounding fluid. This temperature difference is quantified as T p − T f (x p ). The average profile of T p − T f (x p ) zi is computed in vertical slabs of size Δz = 0.004H with central vertical position z i and probability density functions (PDFs) of this quantity are constructed in the BL at the bottom plate (z p < δ T ). From the left panels of fig. 3 we observe that for all values of K the difference between the (average) particle tem-perature and the (average) fluid temperature is indeed increasing with increasing τ T at the bottom plate. The PDFs clearly become wider for larger τ T , again confirming that more extreme temperature differences are found for larger thermal response times as expected. Furthermore, there is an enhanced probability on larger deviations |T p − T f (x p )| for larger values of τ T , when focusing on the left-hand side (lhs) of the PDFs. The temperature difference of the particles with respect to the fluid at x p near the bottom plate is also slightly increasing with increasing K as evident when comparing the top, central and bottom panels of fig. 3. A peak develops on the lhs of the PDFs for increasing K, suggesting that there is indeed a larger probability of larger absolute temperature differences for larger values of K. This is a result of particles with a large thermal expansion coefficient escaping the warm bottom plate region already for a slight temperature increase. Now, only particles that have a much lower temperature with respect to the fluid temperature stay at the plates longer, resulting in a larger absolute temperature difference T p − T f (x p ).

Thermal boundary layer residence time
The thermal response time, τ T , not only influences the temperature difference between particles and the surrounding fluid, but also the time particles reside at the plates before they will escape from the BL due to the buoyancy force. To understand this relation, statistics of the residence time of particles inside the thermal BLs, t δT , are computed for different values of τ T and K, where the resulting PDFs are shown in fig. 4. For τ T 1, the PDFs display a clear peak suggesting that there is a well-defined characteristic time that particles spend inside the thermal BLs. This peak shifts to the right for increasing τ T , so this characteristic residence time increases with increasing τ T as expected. For τ T 1, the PDFs overlap indicating that here t δT is largely independent of τ T . When comparing the different values of K, it is observed that smaller values of t δT are measured for larger values of K, due to particles with a larger thermal expansion coefficient having a quantitatively larger response on temperature fluctuations in the fluid in terms of their mass density.
From the PDFs in fig. 4 it is expected that t δT depends strongly on τ T and K. In fig. 6(a) we show the average residence time of particles inside the thermal BLs at the horizontal plate, t δT , as a function of τ T and for different values of K. Each individual particle can cross the BL multiple times and the average is therefore taken over the total number of times all particles accumulatively cross the BLs. Two regimes can be distinguished in fig. 5, where for τ T 1 the thermal BL residence time is constant, for τ T 1 the values are increasing with increasing τ T . Also, when comparing the three different values of K, the residence time is found to decrease with increasing K. These trends are consistent with the PDFs shown in fig. 4 and confirm that the number of particles inside the BL in fig. 1 is indeed directly related to the time particles spend inside the thermal BL.    Based on the transition from a constant to a ballistic regime, we can estimate the thermal BL residence time as where both a(K) and b(K) are coefficients depending on K. In the limit of small thermal response times, τ T → 0, this equation becomes t δT = a(K) and thus a constant depending only on K. When τ T → ∞ a ballistic behavior t δT = b(K)τ T is found. We perform a fit based on eq. (14) on the DNS data as shown by the dashed lines in fig. 5 and find that the thermal BL residence time indeed depends on τ T as in eq. (14).

Simple 1-dimensional model
To understand in more detail how the dynamics of thermally responsive particles depends on τ T and K, we develop a simple 1-dimensional (1D) model for the thermally responsive particles. The thermal response time, τ T , influences the temperature of particles through the thermal inertia (eq. (5) for the DNS), while the parameter K determines the density ratio β, which is determining the buoyancy force in eq. (4). Therefore, we develop a 1D model for their vertical position in which particles experience both thermal inertia and a buoyancy force: where w is the velocity of particles and we use that T p = T f = T m = 0.5. The velocity, w, is set to zero at z = 0 and z = 1. The fluid temperature T f is now an input parameter of this simple 1D model. In RBC the temperature profile typically shows a large temperature gradient in the thermal BLs, while the temperature in the bulk equals the average temperature [17,36]. Therefore we prescribe a constant mean temperature profile with linear temperature gradients inside the thermal BLs and a constant temperature of T f = 0.5 in the bulk: The thermal BL thickness is set to δ T = 0.022H, equal to the thermal BL thickness measured from the DNS. From the 1D model, we compute the residence time inside the BLs in the same parameter range as in the DNS, 0.05 < τ T < 10 and K = {1.1; 2; 10}, by numerically integrating eqs. (15)-(17) using a second order Adams-Bashforth scheme. Note that the model gives one unique solution and therefore the output is given in terms of t δT and not as an ensemble average t δT as in the DNS. The model results (lines), together with the DNS data (symbols), are shown in fig. 6(a). We observe that the model captures the trend of decreasing t δT with increasing K as observed in the DNS. Also the trend with τ T is recovered where t δT is constant for smaller τ T and is monotonically increasing with τ T for larger τ T .
Let us discuss these two regimes in more detail by looking at the behavior of the model in the limit of i) small thermal response times, τ T → 0, and ii) large thermal response times, τ T → ∞: i) τ T → 0: When the thermal response time is zero, particles are instantaneously adapting their temperature to that of the surrounding fluid such that always T p = T f .
Substituting this into eq. (16) gives us Given that α p > α f and α p T p < 1, we find that in the BL at the top plate where T f > 0.5 the acceleration is negative (a < 0) while in the bottom BL where T f > 0.5 it is positive (a > 0). Equation (18) does not depend on τ T and consequently also the residence time inside the BL (for τ T → 0) is expected to be independent of τ T and to only depend on the thermal expansion coefficient and thus on K. This is exactly what we found in fig. 6(a) for both the 1D model and the DNS in the limit of small τ T . ii) τ T → ∞: For very large thermal response times eq. (17) becomes dT p /dt = 0 and the temperature of particles will be constant and independent of time, such that T p = T p (t = 0). The initial temperature condition is thus fully determining the particle temperature. A particle initially positioned in the bulk will start to move upwards or downwards depending on its initial temperature condition, T p (t = 0). It can be computed that when a particle initially moves upwards, the acceleration and velocity in the top BL are positive such that a particle will end up in the top BL and will get stuck there. Oppositely, an initially downwards moving particle will experience a downwards acceleration and velocity in the bottom BL and will get stuck inside the bottom BL. This means that in the limit of τ T → ∞, t δT → ∞.
These limits are consistent with eq. (14), that was shown to capture the trend of the DNS data. Although not shown here, we verified that the same fitting procedure works for the 1D model confirming that also t δT computed from the 1D model follows eq. (14).
So, both in the model and in the DNS we observe a transition from a constant residence time to a ballistic regime, where the residence time is increasing with increasing thermal response time. However, the transition between the constant and ballistic regimes, occurs at a different value of τ T ≈ 0.1 in the 1D model, compared to the DNS (where the transition occurs around τ T ≈ 1) in fig. 6(a). This might be related to the model being in 1D, while in the DNS particles move in a 3D flow field. As a result the time scales might not be one-to-one comparable. Moreover in the DNS particles are additionally transported towards and away from the plates by the LSC, an effect that is not included in the simple 1D model. Since we do not expect the model and DNS data to match oneto-one, we can re-scale the vertical and horizontal axis for the model results in fig. 6(a). The DNS data, together with the re-scaled data of the 1D model, are shown in fig. 6(b) and now the model matches the DNS results within the error bars. All together we argue that this, though very simple 1D model, captures the trends found in the DNS surprisingly well.

Conclusions
We have studied the dynamics of thermally responsive particles in Rayleigh-Bénard convection. Particles are experiencing both mechanical and thermal inertia, and both fluid and particles exhibit thermal expansion where the thermal expansion coefficient of particles is larger than that of the fluid. Now, particles near the hot bottom plate become lighter than the fluid and particles at the top plate become heavier than the fluid. It is verified that this induces a motion away from the plates, resulting in particles re-suspending from the BLs into the bulk.
This dynamics results in a non-homogeneous distribution of particles throughout the RBC cell. In particular, a regime of thermal response times and thermal expansion coefficients is found where the number of particles at the plate is enhanced compared to a uniform distribution. We have shown that this enhancement is already reached for an increase in the thermal expansion coefficients of particles compared to that of the fluid of ten per cent; K = α p /α f = 1.1. This ratio of thermal expansion coefficients can be achieved in realistic systems, for example gel-like particles in water or oil-water systems.
Upon increasing K, the number of particles at the plates is decreasing, since the particle density responds much stronger to the temperature fluctuations. A regime of large K and small τ T is found where particles escape the BLs almost immediately and where the number of particles inside the thermal BLs is even lower than the uniform distribution. This depletion in the BLs leads to an enhanced number of particles inside the bulk. Increasing τ T has an opposite effect; particles need more time to warm up (cool down) at the bottom (top) plate, increasing the number of particles at the plates and decreasing the number of particles in the bulk.
The number of particles at the plates is expected to depend on the time particles spend inside the thermal BLs at the plates. By quantifying this residence time, t δT , it has been shown that particles do spend a characteristic time inside these BLs, that is moreover depending on τ T and K. In particular, the ensemble average t δT is increasing with decreasing K. For all values of K, t δT is constant for τ T 1 and is increasing with increasing τ T for τ T 1 in the DNS. This trend is confirmed when performing a fit of the function y = a(K) + b(K)x on the DNS data for each value of K. A simple 1D model is developed, where the motion of thermally inertial particles depends exclusively on the buoyancy force, and again both particles and fluid exhibit thermal expansion with α p > α f . This model is shown to capture the trends very well; again the thermal BL residence time is constant for smaller τ T and increasing with increasing τ T for larger τ T , only now the transition occurs at a smaller τ T ≈ 0.1. Also the shift of the curves to lower values of t δT for larger values of K is captured well by the model. When re-scaling the data of the model the DNS and model results match within the error bars, confirming that the model captures the observed trends well. The simple 1D model can thus be used to better understand the interplay between thermal inertia and the buoyancy-driven vertical motion of particles.
We have studied how thermal inertia influences the dynamics of thermally responsive particles, using a pointparticle approach. The dynamics in this point-particle model is already rich and there are many parameters involved. In nature, however, multi-phase fluid systems with different thermal properties for the different phases can become even more complex; e.g. phase transitions in convection in the core of the earth or the presence of deformable vapor bubbles in boiling convection. To study these highly complex systems more advanced numerical techniques, with much higher numerical costs, are necessary. Here we have however shown that DNS with a pointparticle approach is able to give insight into the influence of thermal inertia on the distribution and the temperature statistics of inertial particles in a thermally driven flow where the dispersed phase has different thermal properties than the carrier fluid.