2.1. Particle motion equation

We formulate the governing equation for the particle motion in natural convection over hot surfaces by adding thermophoretic force to the Maxey & Riley (Reference Maxey and Riley1983) equation; the corresponding formulation for forced convection was done by Talbot et al. (Reference Talbot, Cheng, Schefer and Willis1980). For high density, micron sized aerosol particles, as the Basset–Boussinesq history force is much smaller compared with the viscous drag term in the dimensionless Maxey & Riley (Reference Maxey and Riley1983) equation, we neglect this history force (Bergougnoux et al. Reference Bergougnoux, Bouchet, Lopez and Guazzelli2014). Since the particle sizes that we consider are small, Faxen's correction, $(-a^2/6) \nabla ^2u$, to the viscous drag can also be neglected (Maxey & Riley Reference Maxey and Riley1983) in comparison to the terms that are linear in fluid and particle velocities. The governing equation then reduces to

(2.1)\begin{align} \rho_p \frac{{\rm d}\boldsymbol{v}}{{\rm d}t}&=\rho_f\frac{{\rm D}\boldsymbol{u}}{{\rm D}t}-\frac{9}{2}\frac{\mu}{a^2}\frac{1}{C_B}(\boldsymbol{v}-\boldsymbol{u})- \frac{\rho_f}{2}\left(\frac{{\rm d}\boldsymbol{v}}{{\rm d}t}-\frac{{\rm D}\boldsymbol{u}}{{\rm D}t}\right)\nonumber\\ &\quad -\frac{9\rho_f\nu^2}{a^2}C_T\frac{\boldsymbol{\nabla} T}{\langle T_p\rangle}+(\rho_p-\rho_f)\boldsymbol{g}, \end{align}

where the term on the left-hand side is the particle inertia with ${{\rm d}}/{{\rm d}t}= {\partial }/{\partial t} + \boldsymbol{v}\boldsymbol{\cdot } \boldsymbol{\nabla }$ being the material derivative along the particle trajectory and $\boldsymbol{v}$ the particle velocity vector. In the first term on the right in (2.1), ${{\rm D}}/{{\rm D}t} = {\partial }/{\partial t} + \boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }$ is the material derivative along the fluid trajectory, with $\boldsymbol{u}$ being the fluid velocity vector. The second term on the right-hand side of (2.1) is the viscous drag where $C_B$ is any appropriate slip correction factor, such as Basset or Cunningham slip correction factors, to accommodate non-continuum effects, and $\mu$ is the dynamic viscosity of the fluid. The third term on the right-hand side of (2.1) is the added mass, where $\rho _f$ is the fluid density. The next term on the right is the thermophoretic force term, where $T$ is the fluid temperature field and $\langle T_p \rangle$ the undisturbed fluid temperature at the particle centre if the particle was not present at that point. The thermophoretic coefficient in this term is

(2.2)\begin{equation} C_T =\frac{C_s({k_f}/{k_p}+C_tKn)}{(1+3C_mKn)(1+2{k_f}/{k_p}+2C_tKn)}, \end{equation}

with the coefficients $C_m=1.14$, $C_s=1.17$ and $C_t=2.18$ (Talbot et al. Reference Talbot, Cheng, Schefer and Willis1980); $k_f$ and $k_p$ being the fluid and particle thermal conductivities. An alternative expression for $C_T$ in terms of $Kn$ and $k_f/k_p$, given by Beresnev & Chernyak (Reference Beresnev and Chernyak1995), which was shown to be more accurate in the transition regime of Knudsen numbers by Guha & Samanta (Reference Guha and Samanta2014), could also be used instead of (2.2). The last term in (2.1) is the buoyancy acting on the particle, where $\boldsymbol{g}$ is the acceleration vector due to gravity and $\rho _p$ the particle density.

We normalise (2.1) with the characteristic scales near the hot plate in RBC (Theerthan & Arakeri Reference Theerthan and Arakeri1998; Puthenveettil & Arakeri Reference Puthenveettil and Arakeri2005), namely, the near-wall length scale

(2.3)\begin{equation} Z_w = \left(\frac{\nu \alpha}{g\beta \Delta T_w} \right)^{1/3}, \end{equation}

where $\beta$ is the thermal expansivity of the fluid, the near-wall velocity scale

(2.4)\begin{equation} U_w = \frac{\sqrt{\nu \alpha}}{Z_w} \end{equation}

and the near-wall time scale

(2.5)\begin{equation} t_w = Z_w/U_w. \end{equation}

The dimensionless governing equation for particle motion then becomes

(2.6)\begin{equation} \frac{{\rm d}\boldsymbol{v}^*}{{\rm d}t^*}=\frac{3\epsilon}{2+\epsilon}\frac{{\rm D}\boldsymbol{u}^*}{{\rm D}t^*}-\frac{1}{\widehat{St}}(\boldsymbol{v}^*-\boldsymbol{u}^*) - \frac{1}{\widehat{St}} Th \nabla^*{\theta} -\frac{1}{\widehat{St}}Gn\hat{j}, \end{equation}

where all the starred variables denote the dimensionless variables corresponding to the dimensional variables in (2.1), $\epsilon = \rho _f/\rho _p$ and the dimensionless temperature,

(2.7)\begin{equation} \theta = \frac{T-T_\infty}{\Delta T_w}. \end{equation}

Equation (2.6) shows that the following dimensionless numbers determine the particle dynamics. The modified Stokes number, $\widehat {St}$ in (2.6) is

(2.8)\begin{equation} \widehat{St} = \frac{C_B(2+\epsilon)St}{2}, \end{equation}

where the Stokes number

(2.9)\begin{equation} St= \frac{2}{9} \frac{\rho_p a^2/\mu}{t_w} = \frac{2}{9} \frac{1}{\epsilon\sqrt{Pr}} Ra_p^{2/3}, \end{equation}

by using (2.5), (2.3) and (2.4). Here, the particle Rayleigh number,

(2.10)\begin{equation} Ra_p = \left(\frac{g\beta\Delta T_w a^3}{\nu \alpha}\right)^{1/3}, \end{equation}

is the ratio of the time scales of thermal transport by diffusion ($a^2/\sqrt {\nu \alpha })$ to the thermal transport by convection $(\sqrt {\nu \alpha }/g \beta \Delta T_w a)$, both evaluated at the particle length scale. The thermophoretic number,

(2.11)\begin{equation} Th =\frac{1}{Re} \frac{C_TC_B\Delta T_w}{\langle T\rangle}, \end{equation}

in (2.6) represents the dimensionless strength of the thermophoretic effects on the particle, where the Reynolds number

(2.12)\begin{equation} Re = \frac{1}{2} \frac{U_w Z_w}{\nu} = \frac{1}{2\sqrt{Pr}}, \end{equation}

by using (2.3) and (2.4). In (2.6),

(2.13)\begin{equation} Gn = \frac{C_B(1-\epsilon)St}{Fr^2} = \frac{2(1-\epsilon)\widehat{St}}{(2+\epsilon)Fr^2} \end{equation}

is the gravitational number, where the Froude number

(2.14)\begin{equation} Fr = \frac{U_w}{\sqrt{gZ_w}}. \end{equation}

By using (2.8), (2.9) and (2.14) in (2.13),

(2.15)\begin{equation} Gn = \frac{C_B\dfrac{4}{3}{\rm \pi} a^3 (\rho_p-\rho_f)g}{6{\rm \pi} \mu a} \frac{1}{U_w} = \frac{u_{g}}{U_w} = {u^*_{g}}, \end{equation}

where $u^*_{g}$ is the dimensionless gravitational settling velocity and $u_g$ the gravitational settling velocity. Typical values of these dimensionless numbers for olive oil droplets in air are tabulated in table 1 for different diameters of droplets. We note that $\widehat {St} \ll 1$ for the common droplet sizes in aerosols. We have ignored the role of particle interactions and inertial lift forces in the present formulation, effects that could lead to the modification of particle trajectories. Both hydrodynamic and non-hydrodynamic interactions between aerosol particles can play a significant role in the trajectories of aerosol particles, which subsequently influences phenomena like coagulation and deposition (see Patra, Koch & Roy Reference Patra, Koch and Roy2022 and references therein). In the present study, we are interested in the highly dilute limit; thus, the role of interactions is neglected. A particle has a non-zero slip with respect to the local fluid velocity; for small $St$, the deviation is due to the thermophoretic velocity (to be discussed in detail in the next section). Since a particle has a non-zero slip velocity in a local shear flow, it could experience an inertial lift force as was derived in the case of a simple shear flow by Saffman (Reference Saffman1965) (see Candelier, Mehlig & Magnaudet Reference Candelier, Mehlig and Magnaudet2019 for more extensions to more general linear flows). The relative magnitude of the Saffman lift force would depend on the particle Reynolds number defined based on the local shear rate, $Re_{\gamma }$. In the present problem, $Re_{\gamma }=Ra_p^2/\sqrt {Pr}$. From the listed values of $Ra_p$ in table 1, we see $Re_{\gamma }\ll 1$ and thus, the effect of inertial lift is negligible in the present study. This assumption is further corroborated by Yang et al. (Reference Yang, Zhang, Wang, Dong and Zhou2022) where they found the effect of Saffman lift negligible in the study of particle dynamics in Rayleigh–Bénard convection.

Table 1. Values of the various dimensionless numbers for olive oil droplets in air at $Pr=0.7$ for $\Delta T_w \approx 10$ K and $\langle T\rangle \approx 300$ K.

2.2. Limiting cases

When $\widehat {St} \ll 1$, the first term on the right-hand side of (2.6) can be neglected since $1/\widehat {St}$ appears on all the other terms on the right-hand side of (2.6). Then, since $\boldsymbol{v}^* = {\rm d}\kern0.07em\boldsymbol{x}^*/{\rm d}t^*$, for $\widehat {St} \ll 1$, (2.6) reduces to

(2.16)\begin{equation} \widehat{St} \frac{{\rm d}^2\kern0.07em\boldsymbol{x}^*}{{\rm d}{t^*}^2}={-}\left(\frac{{\rm d}\kern0.07em\boldsymbol{x}^*}{{\rm d}t^*}-\boldsymbol{u}^*\right)- Th\nabla^*{\theta} -\widehat{St}\frac{Gn}{\widehat{St}}\hat{j}. \end{equation}

In (2.16), $\widehat {St}$ can approach zero only if $St$ is approaching zero, as shown by (2.8). Hence, when $\widehat {St} \rightarrow 0$, since both $\widehat {St}$ and $Gn$ have a similar dependence on $St$, as evident from (2.8) and (2.13), $Gn/\widehat {St}$ can still remain finite depending on the values that $\epsilon$ and $Fr^2$ take. So when $\widehat {St} \rightarrow 0$ (or $St=0$), the first and the last terms in (2.16) are negligible whereas the second and the third terms remain finite. Further, the term $-Th\nabla ^*{\theta }$ in (2.16) can be re-written by using (2.12), (2.11) and (2.7) as

(2.17)\begin{equation} -Th\nabla^*{\theta}={-}\frac{C_TC_B\Delta T_w}{Re\langle T\rangle}\nabla^*{\theta} = \frac{-2\nu C_T C_B \boldsymbol{\nabla}T}{\langle T\rangle}\frac{1}{U_w} = \frac{\boldsymbol{u}_{t}}{U_w}=\boldsymbol{u}^*_{t}, \end{equation}

where $\boldsymbol{u}_{t}$ is the thermophoretic velocity and $\boldsymbol{u}^*_{t}$ its dimensionless form. Hence, when $\widehat {St} \rightarrow 0$, (2.16) reduces to

(2.18)\begin{equation} \boldsymbol{v}^*=\boldsymbol{u}^* -Th\nabla^*{\theta}=\boldsymbol{u}^* +\boldsymbol{u}^*_{t}, \end{equation}

the zeroth-order equation for the particle velocity. Equation (2.18) implies that when $\widehat {St} \rightarrow 0$, the particle velocity in fluid flows with finite temperature gradients, like thermal boundary layers and plumes, is a vector sum of the fluid velocity and the thermophoretic velocity. Note that this result holds even when $\epsilon \rightarrow 1$ for neutrally buoyant particles, usually employed in particle image velocimetry (PIV) studies. Therefore, even when $St = 0$, neutrally buoyant particles with $\epsilon = 1$ do not follow the fluid trajectory in regions where $\boldsymbol{u}^*_{t}$ is significant. Hence, PIV measurements inside regions like thermal boundary layers and plumes may not be reliable if $\boldsymbol{u}^*_{t}$ values there are not negligible. It is obvious from (2.16) that particles follow the fluid trajectory when $Gn \ll 1$, $Th \ll 1$ and $\widehat {St} = 0$, since $\boldsymbol{v}^* = \boldsymbol{u}^*$ in such a case. In all other cases, the particles would deviate from fluid trajectories under the influence of thermophoretic and gravitational forces, even when $St = 0$.

Since the fluid velocity field and the temperature field that we consider are steady, it can be noted from (2.18) that the particle velocity field also approaches a steady state as $\widehat {St}\rightarrow 0$. This is because in such a case, when $\widehat {St}\rightarrow 0$, the response time required for the particles to reach the velocity of the flow field becomes infinitely small due to the low inertia and hence, the time evolution of the particles in the unsteady state can be neglected in such a case. Then, the initial injection velocities of particles into the flow field become irrelevant in such cases.

2.3. Fluid flow field

We assume that the particles do not modify the flow around them, which is valid for the case of a dilute mixture of small particles, when mass loading of the particle phase is low (Guha Reference Guha2008). It is shown in Appendix A that those conditions are valid in the present study. This assumption allows us to substitute an existing solution of natural convection boundary layers for $\boldsymbol{u}^*(x,y)$ and $\theta (x,y)$ in (2.16) to study the evolution of particles in that flow field. The natural convection boundary layer flow above a horizontal surface has been obtained by Rotem & Claassen (Reference Rotem and Claassen1969), Pera & Gebhart (Reference Pera and Gebhart1973) and Samanta & Guha (Reference Samanta and Guha2012); we here use the similarity solutions proposed by Rotem & Claassen (Reference Rotem and Claassen1969).

Rotem & Claassen (Reference Rotem and Claassen1969) gave the dimensionless stream function,

(2.19)\begin{equation} \psi^*(\eta) = {x^*_{rc}}^{3/5} F(\eta), \end{equation}

and the dimensionless temperature difference,

(2.20)\begin{equation} H(\eta) = \theta = \frac{T-T_\infty}{\Delta T_w}, \end{equation}

as functions of the similarity variable

(2.21)\begin{equation} \eta = \hat{y} {x^*_{rc}}^{{-}2/5}. \end{equation}

Here, the dimensionless $x$ coordinate,

(2.22)\begin{equation} x^*_{rc} = x/L, \end{equation}

where $L$ is the length of the plate. The stretched $y$ coordinate,

(2.23)\begin{equation} \hat{y}=y^*_{rc} Gr_L^{1/5} = \frac{y}{L} Gr_L^{1/5}, \end{equation}

where the Grashoff number

(2.24)\begin{equation} Gr_L = \frac{g\beta \Delta T_w L^3}{\nu^2} = \frac{Ra_L}{Pr}, \end{equation}

where $Ra_L = (-g\beta \Delta TL^3)/\nu \alpha$ is the Rayleigh number based on $L$, so that, using (2.21), (2.22), (2.23) and (2.24),

(2.25)\begin{equation} \eta = \frac{y}{L}\left(\frac{x}{L}\right)^{{-}2/5}\left(\frac{Ra_L}{Pr}\right)^{1/5}= \frac{y}{(x^2Z_w^3)^{1/5}}\frac{1}{Pr}. \end{equation}

From Rotem & Claassen's (Reference Rotem and Claassen1969) stretched horizontal and vertical velocity components, given in (D3) and (D4), we obtain the dimensionless velocity components for the present study as

(2.26)\begin{equation} u_x^* = \frac{\partial \psi^*}{\partial \hat{y}} = Pr^{1/10} {x^*}^{1/5} F'(\eta) \end{equation}

and

(2.27)\begin{equation} u_y^* ={-} \frac{\partial \psi^*}{\partial {x_{rc}}^*} = \frac{1}{5} Pr^{3/10} {x^*}^{{-}2/5} (2\eta F'(\eta) - 3F(\eta)), \end{equation}

by using (D1), (B1), (B2), (2.3) and (2.4).

Values of $F(\eta )$ from (2.19) and $H(\eta )$ from (2.20) for various $Pr$ are given by Rotem & Claassen (Reference Rotem and Claassen1969) by numerical solution of the natural convection boundary layer equations (2.28). The values of $\eta$, and hence the values of the velocity components (2.26), (2.27) and the temperatures (2.20) at any location inside the boundary layer are then determined by the values of $Ra_L$ and $Pr$ through (2.25). Here, $F(\eta )$ and $H(\eta )$ are obtained by numerically solving the following boundary layer equations of Rotem & Claassen (Reference Rotem and Claassen1969),

(2.28a)$$\begin{gather} 5F'''+3FF''-{F'}^2 = 2G - 2\eta G', \end{gather}$$

(2.28b)$$\begin{gather}H = G' \end{gather}$$

and

(2.28c)\begin{equation} H'' + \tfrac{3}{5}Pr FH'=0, \end{equation}

where $'$ indicates derivative with respect to $\eta$. Additionally,

(2.29)\begin{equation} G(\eta) = \left(\frac{p-p_{\infty}}{\rho_{\infty} \nu^2} L^2 + \frac{g{y^*_{rc}}L^3}{\nu^2}\right) Gr^{{-}4/5} {x^*_{rc}}^{{-}2/5} \end{equation}

is the dimensionless pressure function, where $p$ is the local fluid pressure, $p_{\infty }$ the pressure of the ambient fluid far away from the plate and $\rho _{\infty }$ the corresponding density. The boundary conditions for solving (2.28) include

(2.30a)\begin{equation} F=F'=0,\quad H=G'=1 \quad \mathrm{at} \ \eta =0 \end{equation}

and

(2.30b)\begin{equation} F'=H=G=0 \quad \mathrm{when} \ \eta \rightarrow \infty. \end{equation}

2.4. Numerical method

Numerical solutions of (2.16) were obtained for 5000 particles, initially placed randomly in the domain with zero initial velocity at $\widehat {St} = 10^{-5}$ for different combinations of $Gn$ and $Th$. The horizontal domain started at $x = 0.01Z_w$ (or $x^* = 0.01$) to exclude the origin of the boundary layer, which is a singular point. The domain extended horizontally up to $50Z_w$, approximately the average length of a growing boundary layer before it turns into a plume (Theerthan & Arakeri Reference Theerthan and Arakeri1998). In the vertical direction, since the boundary layer thickness $\delta \sim Z_w$, the domain extended from $y=0$ to $y=40Z_w$ to ensure that both the thermal and the velocity boundary layers were well within the domain.

To maintain a constant number of particles in the domain, whenever a particle exited, another particle was introduced through the top of the domain at a random $x$ position with zero velocity. As discussed in § 2.2, the initial velocities of the particles become insignificant when $\widehat {St} \ll 1$, as they reach the surrounding fluid velocity in a negligible amount of time. Therefore, the initial velocities at the start of the computation and the injection velocities of particles through the top boundary were set to zero for convenience.

Equation (2.16) was numerically solved using the fifth-order Runge–Kutta method with adaptive step size selection, implemented via the ODE45 tool in MATLAB$\circledR$. The particles evolved in the steady flow field described by (2.26), (2.27) and (2.20). To solve (2.16) for particle positions, $u_x^*$ in (2.26), $u_y^*$ in (2.27) and $\theta$ in (2.20) were obtained by interpolating the solutions of the boundary layer equations (2.28). For efficient implementation of interpolation, (2.16) was integrated using a curve-fitted form of $F(\eta )$ and $H(\eta )$, which were used to obtain $u^*_x$ in (2.26), $u^*_y$ in (2.27) and $\theta$ in (2.20). The fits for $F(\eta )$ and $H(\eta )$ are shown in figure 2.

Figure 2. Comparison of the curve fits obtained for $F(\eta )$ and $H(\eta )$ to the data points obtained from Rotem & Claassen (Reference Rotem and Claassen1969) for $Pr=0.7$. Data points from Rotem & Claassen (Reference Rotem and Claassen1969) for: $\triangledown$, $F(\eta )$; and red $\triangledown$, $H(\eta )$. Curve fits for: black – – –, $F(\eta )$; and blue – – –, $H(\eta )$.

The ODE45 algorithm automatically selects smaller time steps where the solutions change abruptly, ensuring the accuracy of the solution based on the specified tolerance. The current simulations were conducted with a relative tolerance of $10^{-6}$. Solving the trajectory of a particle from one end of the domain to the other took approximately 5 minutes. Since each particle trajectory could be independently solved using multiple threads in parallel, the computational time was drastically reduced.

The curve fits in figure 2 for $Pr = 0.7$ were used only to obtain the particle dynamics discussed in § 3. The scaling laws in the subsequent sections were derived using $F(\eta )$ and $H(\eta )$ obtained from numerical solutions of the boundary layer equations in (2.28).

3.1. Flow streamlines and particle trajectories

We analysed the dynamics of particles inside the boundary layer at $\widehat {St}=10^{-5}$ for different combinations of $Gn$ and $Th$. Figure 3(a) shows the particle positions and fluid's streamlines, while figure 3(b) shows the particle trajectories and fluid streamlines corresponding to $Gn=0$ and $Th=1$ at $t^*=45$. Figure 3(c,d) shows the same at ${Gn=0.14}$ and $Th=1$ at $t^*=60$. The corresponding motion of particles can be clearly viewed in the supplementary movies for the cases with and without gravity. The particle trajectories in figure 3(b,d) and the motion of the particles in the two videos corresponding to them show that the particles at the top of the domain are initially dragged down by the entrainment fluid flow into the boundary layer. For the low-gravitational-force situations, as in figure 3(a), once the particles enter the boundary layer, they mostly take a right turn before finally exiting the domain through the right end. For cases when $Gn$ is non-negligible compared with $Th$, as in figure 3(c), many of the particles settle on the plate, as they move to the right; the location of settling depends on the relative strength of the gravitational force that the particles experience.

Figure 3. For $Gn=0$ and $Th=1$: (a) instantaneous particle positions and fluid streamlines; (b) fluid streamlines and particle trajectories. (c,d) The same for $Gn=0.14$ and $Th=1$; blue ——, fluid stream lines; blue – – – –, velocity boundary layer edge; – – – –, particle trajectories; Please refer to the particle motion videos corresponding to panel (a) at $Gn=0$ and to panel (c) at $Gn=0.14$, provided in the supplementary movies, to see the evolution of particle positions.

In figure 3(a,c), we see a dust-free region within the boundary layer, which forms under the influence of thermophoretic force. Both the figures show the final steady-state situation, after the initial transient development of the dust-free region is completed, which takes a dimensionless time of approximately $t^*\sim 50$. In figure 3(a), for the $Gn=0$ case, we see a dust-free region that keeps growing in height as we move along the positive $x$ direction; there is no settling of particles on the hot surface. For this case, the height of the dust-free region at any given $x$ depends on the particular value of $Th$. In contrast, the height of the dust-free region in figure 3(c), after an initial increase, keeps decreasing as $x$ increases. In the finite $Gn$ case, the dust-free region ends at a finite distance from the origin of the boundary layer, beyond which the particles start to settle under the influence of gravity; the length of this dust-free region in figure 3(c) is approximately $x^* = 45$. We observed that if $Gn$ is increased, the particles start to settle at a shorter $x^*$ making the dust-free region smaller in size and height; increasing $Th$ was seen to have an opposite effect. In both the figures, particles that reach a region close to the leading edge of the boundary layer turn towards the left and exit the domain through the left end of the domain.

Since the particle flow divides into left and right near the leading edge, there exists a saddle point close to the leading edge in the flow field, characterised by zero particle velocity and intersecting particle trajectories.

We obtained such a saddle point in the flow field by solving (2.16) numerically and obtaining the co-ordinates where the particle velocity goes to zero; the intersecting particle trajectories were then obtained by solving (2.16) either backwards in time or forwards in time starting from any point infinitesimally close to the saddle point. Figure 4 shows the obtained particle trajectories close to the saddle point for the case in figure 3(a), where the dashed lines show the trajectories of the particles and the solid lines showing the separatrices that pass through the saddle point. At the saddle point, four particle trajectories, two of which end there starting from the top and the bottom of the domain, and two leaving the saddle point towards the left and right so as to exit the domain can be observed in figure 4; we call these trajectories separatrices. For the $Gn=0$ case shown in figure 4, the separatrix leaving to the right grows in height and forms the boundary of the dust-free region. For the corresponding $Gn=0.14$ case, as we show later in figure 5, the height of the separatrix leaving to the right would diminish as $x$ increases and the separatrix would eventually touch the plate surface at a finite distance from the origin, making the dust-free region finite in size. For such cases, all trajectories of the particles meet the plate at finite distances from the origin.

Figure 4. Particle trajectories close to the saddle point for the case shown in figure 3(a) with $Gn=0$ and $Th=1$. The solid lines show the four separatrices with arrows showing the direction of the particle motion on the separatrix. The dashed lines show the trajectories of the other particles in the flow field. Particle positions in the domain at a later instant after the dust-free region has formed are shown with black dots.

Figure 5. Separatrices for varying ratios of $Gn$ to $Th$ for $Th=1$ and $Th=10$ for $\widehat {St}=10^{-5}$. The blue set of curves corresponds to $Th=1$ and the black set of curves corresponds to $Th=10$. ——, $Gn/Th=0$; – – – –, $Gn/Th=0.001$; $\cdots \cdots$, $Gn/Th=0.1$; and -$\cdot$-$\cdot$-$\cdot$-, $Gn/Th=0.14$. The red dashed lines correspond to $y^* \sim {x^*}^{2/5}$, the scaling of natural convection boundary layer thickness.

The separatrices then divide the flow into four sub regions, labelled as I to IV in figure 4, which are characterised by a different behaviour of particles. In region I, the particles coming down from the top take a right turn and move in the positive $x$ direction. In region II, particles coming down take a left turn and move in the negative $x$ direction. In regions III and IV, the particles close to the plate are lifted off due to the upward component of the thermophoretic force and then take a left or a right turn, respectively. As is clear from the trajectories shown in figure 4, all the particles in regions III and IV either exit from the domain through the right/left or settle on the plate; there is no supply of particles to these regions from the regions I and II on the top. The regions III and IV will hence evolve with time into a dust-free region, with the separatrix AOB as its top boundary; the separatrix AOB is hence the steady-state boundary of the dust-free region close to the hot surface.

3.2. Particle velocity for small and finite $\widehat {St}$

For finite but small $\widehat {St}$, for which case the present simulations are conducted, the above flow pattern with saddle point and separatrices can be better understood based on the simplification of (2.16). For $\widehat {St} \ll 1$, but not equal to zero, (2.18) can be modified as

(3.1)\begin{equation} \boldsymbol{v}^*=\boldsymbol{u}^* +\boldsymbol{u}^*_{t}+\widehat{St}\boldsymbol{a^*}, \end{equation}

where $\boldsymbol{a}^*$ is an inertial correction that needs to be evaluated. Substituting (3.1) in (2.6), neglecting all terms of order higher than $\widehat {St}$ and then using (2.15), we obtain the inertia corrected particle velocity as

(3.2)\begin{equation} \boldsymbol{v}^*=\boldsymbol{u}^* +\boldsymbol{u}^*_{t}+\boldsymbol{u}^*_{g}-\widehat{St} \left(\frac{2(1-\epsilon)}{2+\epsilon}\boldsymbol{u}^*\boldsymbol{\cdot }\boldsymbol{\nabla}{\boldsymbol{u}^*} +\boldsymbol{u}^*\boldsymbol{\cdot }\boldsymbol{\nabla}{\boldsymbol{u}^*_{t}}+\boldsymbol{u}^*_{t}\boldsymbol{\cdot }\boldsymbol{\nabla}{\boldsymbol{u}^*}+\boldsymbol{u}^*_{t}\boldsymbol{\cdot }\boldsymbol{\nabla}{\boldsymbol{u}^*_{t}}\right). \end{equation}

Equation (3.2) is a modification of the slow manifold (Sapsis & Haller Reference Sapsis and Haller2008) of Maxey & Riley (Reference Maxey and Riley1983) equations with thermophoretic effects. When $\widehat {St} \ll 1$, the last term in (3.2) is negligible while the third term on the right-hand side of (3.2) may still be significant. Hence, for $\widehat {St} \rightarrow 0$ but finite, (3.2) simplifies to

(3.3)\begin{equation} \boldsymbol{v}^* = \boldsymbol{u}^* +\boldsymbol{u}^*_{t} +\boldsymbol{u}^*_g, \end{equation}

implying that the particle velocity will be a superposition of the fluid velocity, the thermophoretic velocity (2.17) and the Stokes settling velocity (2.15).

The presence of a saddle point and separatrices in the flow field can now be better understood in terms of (3.3). Since $\boldsymbol{u}_{t}^*$ acts opposite to the direction of the temperature gradient, the horizontal and the vertical components of $\boldsymbol{u}_{t}^*$, namely ${u_{t_x}^*}$ and $u_{t_y}^*$, expressions for which are derived in Appendix B, will be in the negative $x$ direction and in the positive $y$ direction, respectively. Since the entraining fluid coming down towards the plate takes a turn in the positive $x$ direction inside the boundary layer, the horizontal component $u^*_x$ of the fluid velocity in the boundary layer will be in the positive $x$ direction; $u_y^*$ will be in the negative $y$ direction. Here, $\boldsymbol{u}_g^*$ has only a vertical component in the negative $y$ direction in the case of a horizontal plate. The saddle point is the point where $u_{t_y}^*$ balances $u_y^*$ and $u_g^*$, as well as where $u_{t_x}^*$ balances $u^*_x$. Then, in regions I and IV, $u^*_x$ is more than $u_{t_x}^*$, and similarly, in regions II and III, $u_{t_x}^*$ is more than $u^*_x$. The thermophoretic velocity $\boldsymbol{u}^*_t$ hence balances the components of $\boldsymbol{u}^*$ and $\boldsymbol{u}^*_g$ to create a saddle point, with the four separatrices in the flow field separating regions where the different forces dominate.

Figure 5 shows the plots of separatrices for four different values of $Gn/Th$ at $Th=1$ and $Th=10$. As expected, for both values of $Th$, the areas under the separatrix, which are the dust-free regions, reduce as $Gn$ increases and the gravitational force becomes more dominant. The dust-free region shrinks to a small patch close to the origin of the boundary layer when $Gn$ is of the same order as $Th$. Figure 5 shows that when $Gn=0$ or when $Gn \ll Th$, the height of the dust-free region, at distances away from the saddle point, scales as $y^* \sim {x^*}^{2/5}$, shown as the red dashed line, in the same manner as the natural convection boundary layer thickness. However, as is obvious from figure 3(a), the thickness of the dust-free region is much lower than the boundary layer thickness. Further, the thickness of the dust-free region is also a function of $Th$, as shown in figure 5, unlike the boundary layer thickness. As we saw earlier, since the separatrix AOB is the boundary of the dust-free region, an expression obtained for the separatrix would also be an equation for the dust-free region that would expectedly capture the dependence of the dust-free region on the fluid and the flow properties. We now obtain such an expression for the separatrix when the Stokes settling velocity $\boldsymbol{u}^*_g \ll \boldsymbol{u}^*_{t}$, the thermophoretic velocity, which occurs when $Gn \ll Th$.

The present novel method of finding the dust-free region boundary, by finding the equation for the separatrix, was chosen as earlier proposed conditions to obtain dust-free regions failed in the present situation. Stratmann et al. (Reference Stratmann, Fissan, Papperger and Friedlander1988) used the locus of all points having zero vertical velocity as the dust-free region boundary; however, such a curve was observed to be a constant $\eta$ curve, different from the dust-free region curve, in the present case. The locus of the maximum negative values of the divergence of the particle velocity field (3.3) was also considered as a criterion to identify the boundary of the dust-free region. This criterion was considered since the particles move out of the dust-free region to the boundary of the dust-free region, and due to its similarity with the recently proposed condition by Shevkar et al. (Reference Shevkar, Vishnu, Mohanan, Koothur, Mathur and Puthenveettil2022) to identify plumes, which are converging flow regions. However, the locus of all such points was found to be another constant $\eta$ line, which did not coincide with the dust-free region boundary.

3.3. Scaling of the dust-free region height when $Gn \ll Th$ and $\widehat {St}\ll 1$

As we saw in § 3.2, when $Gn \ll Th$, the dust-free region height has a similar dependence on $x^*$ as the boundary layer thickness, but which, in contrast to the boundary layer thickness, depends on $Th$. For such cases, an instance of which is shown in figure 4, we now find an expression for the separatrix AOB; this gives the scaling of the dust-free region height for $Gn \ll Th$. For $\widehat {St} \ll 1$, by taking the ratio of $y$ and $x$ components of (3.3), we obtain a differential equation for the general trajectory of a particle as

(3.4)\begin{equation} \frac{{{\rm d}\kern0.05em y}^*}{{{\rm d}\kern0.07em x}^*}=\frac{v_y^*}{v_x^*}=\frac{u_y^* +u_{t_y}^*+u_{g}^*}{u_x^* +u_{t_x}^*}. \end{equation}

From (2.15) and (2.17), $u_{g}^*$ in (3.4) can be neglected in comparison to $u_{t_y}^*$ since $Gn \ll Th$, to obtain

(3.5)\begin{equation} \frac{{{\rm d}\kern0.05em y}^*}{{{\rm d}\kern0.07em x}^*} \approx \frac{u_y^* +u_{t_y}^*}{u_x^* +u_{t_x}^*}. \end{equation}

As we saw in § 3, $u^*_{t_x}$ is expected to be smaller than $u^*_x$ in regions I and IV, and hence, also along the separatrix OB. The expressions for $u^*_x$ and $u^*_{t_x}$ given by (2.26) and (B7) further show that for a given $Th$ on constant $\eta$ line, like the separatrix, $u^*_{t_x}/u^*_x \sim {x^*}^{- 6/5}$. Hence, for a given $Th$ far away from O, since the separatrix follows an $\eta =$ constant variation along OB in figure 4, $u^*_{t_x}$ can be neglected in comparison to $u^*_x$; (3.5) then reduces to

(3.6)\begin{equation} \frac{{{\rm d}\kern0.05em y}^*}{{{\rm d}\kern0.07em x}^*}=\frac{v_y^*}{v_x^*}=\frac{u_y^* +u_{t_y}^*}{u_x^*}. \end{equation}

Substituting (2.26), (2.27) and (B8) for $u_x^*$, $u_y^*$ and $u_{t_y}^*$ in (3.6) results in

(3.7)\begin{equation} \frac{{{\rm d}\kern0.05em y}^*}{{{\rm d}\kern0.07em x}^*}=\frac{({-}3F(\eta)+2\eta F'(\eta)-5ThH'(\eta)Pr^{{-1}/{2}})Pr^{{1}/{5}}}{5F'(\eta){x^*}^{{{3}/{5}}}}. \end{equation}

We now integrate (3.7) with respect to $x^*$ for $\eta =$ constant, a property of the separatrix for $Gn \ll Th$ (see figure 5), to obtain the equation for the separatrix at any $Th$, when $Gn \ll Th$ and $\widehat {St} \ll 1$, as

(3.8)\begin{equation} y^*=\frac{({-}3F(\eta)+2\eta F'(\eta)-5Pr^{-{1}/{2}}ThH'(\eta)){x^*}^{{2}/{5}}Pr^{{1}/{5}}}{2F'(\eta)}. \end{equation}

Substituting $\eta = {x^*}^{-{2}/{5}} y^* Pr^{-{1}/{5}}$ from (B2), (3.8) simplifies to

(3.9)\begin{equation} F(\eta)=\frac{-5}{3\sqrt{Pr}}ThH'(\eta), \end{equation}

an equation for the separatrix at any $Th$ and $Pr$ when $Gn \ll Th$ and $\widehat {St}\ll 1$. Using the expression for $Re$ from (2.12) in (2.11), and defining a modified thermophoretic number,

(3.10)\begin{equation} Th_m = \frac{Th}{\sqrt{Pr}} = \frac{2C_T C_B \Delta T_w}{\langle T\rangle}, \end{equation}

as a $Pr$ independent dimensionless number for the thermophoretic force, (3.9) can be rewritten as

(3.11)\begin{equation} F(\eta)={-}\tfrac{5}{3}Th_mH'(\eta). \end{equation}

Since $F(\eta )$ and $H'(\eta )$ are functions of $Pr$, the solution of (3.11) will be a function of $Th_m$ and $Pr$. Using the values of $F(\eta )$ and $H'(\eta )$ obtained from the numerical solution of (2.28a) to (2.28c) for any $Pr$, in (3.11), for any given value of $Th_m$, we get the value of $\eta$ which specifies the separatrix OB in figure 4, at distances far from the saddle point O. We will hereinafter refer to these values as $\eta _{df}$, the subscript $df$ indicating that these curves specify the boundary of the dust-free region. Figure 6 shows the separatrices obtained from the numerical solution of (2.16) for $Th=1$, $Th=10$ and $Th=100$ at $Pr=0.7$, $Gn=0$ and $\widehat {St}=10^{-5}$, along with the $\eta _{df}$ values predicted by (3.11) for the same values of $Th$ and $Pr$. It is clear from the figure that the curves of separatrices predicted by (3.11) approximate the dust-free region fairly well in the entire domain, and more accurately, in the regions away from the saddle point. The plots of $\eta _{df}$ versus $Th_m$ for $Pr$ ranging from $0.1$ to $10$, obtained from the solution of (3.11), are shown in figure 7. It is observed from the figure that for smaller $\eta _{df}$ ($\eta _{df}\le 1$), $\eta _{df}$ increases as $Pr$ increases while for larger $\eta _{df}$ ($\eta _{df}\ge 3$), $\eta _{df}$ decreases as $Pr$ increases. Figure 7 also shows that $\eta _{df}$ has different scaling dependences on $Th_m$ for small and large $\eta _{df}$; we now obtain these scaling laws for small and large $\eta$, with the $Pr$ dependence obtained later in § 4.

Figure 6. Comparison of separatrices from the numerical solution of (2.16) with the solution of (3.11) for $\widehat {St} = 10^{-5}$, $Gn=0$ and $Pr=0.7$, for different $Th$. Numerical solution of (2.16) for: -$\cdot$-$\cdot$-, $Th=1$; ——, $Th=10$; – – –, $Th=100$. The solution of (3.11) for: –$\diamond$–, $Th=1$; –$\triangle$–, $Th=10$; –$\square$–, $Th=100$.

Figure 7. Variation of dimensionless dust-free region height $\eta _{df}$ with the modified thermophoretic number (3.10) for various $Pr$: red $\circ$, $Pr=0.1$; $*$, $Pr=0.3$; orange $\diamond$, $Pr=0.5$; $\triangledown$, $Pr=0.7$; green star, $Pr=1$; red $\triangleright$, $Pr=3$; $+$, $Pr=7$; blue $\square$, $Pr=5$; brown $\triangle$, $Pr=10$; – – –, $\eta _{df} \sim \sqrt {Th_m}$; $\cdots \cdots \cdots$, $\eta _{df} \sim \ln (Th_m)$. The inset shows the curves $-(5/3)Th_mH'(\eta )$ for; blue – – –, $Th_m =1$; green – – –, $Th_m = 5$; magenta – – –, $Th_m = 10$; and ——, $F(\eta )$ for $Pr=0.7$.

3.3.1. Limit $\eta \rightarrow 0$

We obtain the dependence of $\eta _{df}$ on $Th_m$ for smaller values of $\eta _{df}$ by using approximations of $F'(\eta )$ and $H(\eta )$ in (3.11) as $\eta \rightarrow 0$. A local analysis of (2.28a) to (2.28c) around $\eta =0$, given in Appendix C.1, shows that $F'(\eta )$ and $H(\eta )$ are linear functions in $\eta$ as $\eta \rightarrow 0$ (see (C8)). We hence approximate $F'(\eta ) = 2 M(Pr) \eta$ and $H(\eta ) = 1- N(Pr) \eta$ so that

(3.12a)\begin{equation} F(\eta) = M(Pr) \eta^2 \end{equation}

and

(3.12b)\begin{equation} H'(\eta) ={-}N(Pr), \end{equation}

where $M(Pr)$ and $N(Pr)$ are some positive real functions of $Pr$. Using (3.12a) and (3.12b) in (3.11), we obtain

(3.13)\begin{equation} \eta_{df} = Q(Pr) \sqrt{Th_m}, \end{equation}

where using (2.22), (2.23) and (2.3) in (2.21),

(3.14)\begin{equation} \eta_{df} = \frac{y_{df}}{(Z_w^3 x^2 Pr)^{1/5}} \end{equation}

is the dimensionless dust-free region height, with $y_{df}$ being the dimensional dust-free region height and

(3.15)\begin{equation} Q(Pr) = \sqrt{\frac{5}{3} \frac{N(Pr)}{M(Pr)}} \end{equation}

is a function of $N(Pr)$ and $M(Pr)$. Since both $N(Pr)$ and $M(Pr)$ are positive real functions, $Q(Pr)$ is always positive for a hot plate, because of which, for all $Pr$, $\eta _{df}$ always increase as $\sqrt {Th_m}$ for small $\eta$, as shown in figure 7.

3.3.2. Limit $\eta \rightarrow \infty$

The dependence of $\eta _{df}$ on $Th_m$ for large values of $\eta _{df}$ is obtained as follows. As $\eta \rightarrow \infty$, the boundary condition (2.30b) implies that

(3.16a)\begin{equation} F(\eta) = C_1(Pr), \end{equation}

where $C_1(Pr)$ is a positive real function of $Pr$. Using (3.16a) in (2.28c) and then integrating gives

(3.16b)\begin{equation} H'(\eta)= C_2(Pr)\,{\rm e}^{{-}3C_1(Pr)Pr\eta/5}, \end{equation}

where $C_2(Pr)$ is a negative real function of $Pr$. Using (3.16a) and (3.16b) in (3.11) gives

(3.17)\begin{equation} \eta_{df} = R(Pr) + S(Pr) \ln(Th_m), \end{equation}

where

(3.18a,b)\begin{equation} R(Pr)= \frac{-5}{3PrC_1(Pr)}\ln\left(\frac{-3C_1(Pr)}{5C_2(Pr)}\right) \quad \mathrm{and} \quad S(Pr) = \frac{5}{3} \frac{1}{C_1(Pr) Pr} \end{equation}

are functions of $Pr$. Since $S(Pr)$ is a positive real function in $Pr$, we see that $\eta _{df}$ increases as $\ln (Th_m)$ for all $Pr$ at large $\eta$ as shown in figure 7; a higher $\eta _{df}$ always corresponds to a higher value of $Th_m$ for a given $Pr$. This dependence could also be explained from a graphical depiction of (3.11). The inset of figure 7 shows a plot of $-(5/3)Th_mH'(\eta )$, the right-hand side of (3.11), versus $\eta$ for $Pr=0.7$ along with the plot for $F(\eta )$, the left-hand side of (3.11), versus $\eta$ for the same $Pr$. The curves for $-(5/3)Th_mH'(\eta )$ corresponding to $Th_m=1$, $Th_m=5$ and $Th_m=10$ cross the curve for $F(\eta )$ at higher values of $\eta$, as $Th_m$ is increased from 1 to 10. In (3.13) and (3.17), the $Pr$ dependence of $\eta _{df}$ is not fully specified; we now obtain the complete expressions for $\eta _{df}$ in terms of $Pr$ and $Th_m$ for different ranges of $Pr$.

4.1. Scaling for $Pr\gg 1$

4.1.1. Scaling for small $\eta$

Rotem & Claassen (Reference Rotem and Claassen1969) provides the inner boundary layer equations for the case of $Pr \gg 1$ as follows:

(4.1a–c)\begin{equation} 5\tilde{F}'''_{1}= 2\tilde{G}_{1} - 2\tilde{\eta} \tilde{G}'_{1}, \quad \tilde{H}_{1} = \tilde{G}'_{1}, \quad \tilde{H}''_{1} + \tfrac{3}{5} \tilde{F}_{1}\tilde{H}'_{1}=0, \end{equation}

where ${\tilde {F}}_{1}$, ${\tilde {G}}_{1}$ and ${\tilde {H}}_{1}$ are the dimensionless velocity, pressure and temperature functions obtained by using appropriate transformations as in (D6) and (D7). The boundary conditions include

(4.2a)\begin{equation} \tilde{F}_{1}=\tilde{F}'_{1}=0; \quad \tilde{H}_{1}=\tilde{G}'_{1}=1 \quad \mathrm{at} \ \tilde{\eta}_{1} =0 \end{equation}

and

(4.2b)\begin{equation} \tilde{F}''_{1}=0; \quad \tilde{H}_{1}=\tilde{G}_{1}=0 \quad \mathrm{when}\ \tilde{\eta}_{1} \rightarrow \infty, \end{equation}

where $\tilde {\eta }_{1}$ is the similarity variable defined for $Pr\gg 1$ in (D6). Since there is no $Pr$ appearing in (4.1), $\tilde {F}_{1}$ and $\tilde {H}_{1}$ will only depend of $\tilde {\eta }_{1}$ and not on $Pr$.

A local analysis of (4.1) near $\tilde {\eta }_{1}=0$, given in Appendix C.2, shows that $\tilde {F}'_{1}(\tilde {\eta }_{1})$ and $\tilde {H}_{1}(\tilde {\eta }_{1})$ are linear when $\tilde {\eta }_{1} \rightarrow 0$. In such a case,

(4.3)\begin{equation} \tilde{F}_{1}(\tilde{\eta}_{1})=C_1\tilde{\eta}^2_{1}, \end{equation}

where $C_1=\tilde {F}''_{1}(\tilde {\eta }_{1}=0)/2=0.49$ and

(4.4)\begin{equation} \tilde{H}'_{1}(\tilde{\eta}_{1})=C_2, \end{equation}

where $C_2=\tilde {H}'_{1}(0)=-0.46$. It is to be noted that $\tilde {F}_{1}(\tilde {\eta }_{1})$ in (4.3) and $\tilde {H}'_{1}(\tilde {\eta }_{1})$ in (4.4) are independent of $Pr$, since the boundary layer equations (4.1) do not have $Pr$ in them. So $C_1$ and $C_2$ are just constants and not functions of $Pr$. From (D9), (D3) and (D5b), we have

(4.5)\begin{equation} F'(\eta) = Pr^{{-}3/5}\tilde{F}'_{1}(\tilde{\eta}_{1}). \end{equation}

Similarly, from (D5a), (D6c) and (2.21), we obtain

(4.6)\begin{equation} \eta=Pr^{{-}1/5}\tilde{\eta}_{1}. \end{equation}

Using (4.6) in (4.5) and integrating, since $\tilde {F}_{1}(0)$ and $F(0)$ are equal to zero gives us

(4.7)\begin{equation} F(\eta) = Pr^{{-}4/5}\tilde{F}_{1}(\tilde{\eta}_{1}). \end{equation}

Differentiating (D6b) and using (4.6), we get

(4.8)\begin{equation} H'(\eta)=Pr^{1/5}\tilde{H}'_{1}(\tilde{\eta}_{1}). \end{equation}

Even though both $\tilde {F}_{1}(\tilde {\eta }_{1})$ and $\tilde {H}'_1(\tilde {\eta }_{1})$ are independent of $Pr$ (see (4.3) and (4.4)), the transformation back to the original boundary layer variables $F(\eta )$ and $H'(\eta )$ in (4.7) and (4.8) brings back the $Pr$ dependence of the velocity and the temperature fields. Using (4.7) and (4.8) to convert (4.3) and (4.4) to be in terms of $\eta$, then substituting in (3.11) and using (4.6), we get

(4.9)\begin{equation} \eta_{df}=1.25 \sqrt{Th_m} Pr^{3/10}. \end{equation}

Equation (4.9) is the theoretically obtained scaling relation for $\eta _{df}$ in the small $\eta$ limit for $Pr \gg 1$. A comparison of (3.13) with (4.9) makes it clear that the unspecified $Pr$ dependence in (3.13) is

(4.10)\begin{equation} Q(Pr)=1.25 Pr^{3/10}. \end{equation}

In figure 8, we compare (4.9) with $\eta _{df}$ obtained from (3.11) for $Pr=2$, $Pr=10$ and $Pr=100$; where $F(\eta )$ and $H'(\eta )$ corresponding to $Pr=2$, $Pr=10$ and $Pr=100$ used in (3.11) were obtained by numerically solving (2.28). Since (4.9) is valid for $Pr \gg 1$, we expect the data of $Pr=100$ to approach (4.9) in the small $\eta$ limit, while the data of $Pr=2$ to deviate. Surprisingly, in figure 8, the data of $Pr=2$, $Pr=10$ and $Pr=100$ collapse onto (4.9) suggesting that (4.9) is valid not just for $Pr \gg 1$, but also for any $Pr>1$, for small $\eta$.

Figure 8. Scaling of the dimensionless dust-free region height $\eta _{df}$ for $Pr \gg 1$ in the small and the large $\eta$ limits. Numerical solution of (3.11) for: blue $\circ$, $Pr=2$; green $\triangle$, $Pr=10$; and red $\diamond$, $Pr=100$. The scaling relations ———, (4.9) and – – – –, (4.13). The inset shows the variation of $\tilde {F}(\tilde {\eta })$ with $\tilde {\eta }$ for $Pr \gg 1$, obtained by solving the inner boundary layer equations (4.1).

4.1.2. Scaling for large $\eta$

For $Pr \gg 1$, the thermal boundary layer is an inner boundary layer; there is negligible temperature gradient above it in the velocity boundary layer. Since thermophoresis and the resultant formation of the dust-free region are due to temperature gradients, for $Pr \gg 1$, we expect the dust-free region to not extend beyond this inner thermal boundary layer. The scaling of the dust-free region height in the large $\eta$ limit has to be then obtained from an analysis of the top regions of the inner thermal boundary layer. We hence look at the inner boundary layer equations (4.1), with the corresponding boundary conditions (4.2) to find the large $\eta$ scaling of $\eta _{df}$ at $Pr\gg 1$.

From the $\tilde {\eta }_{1} \rightarrow \infty$ boundary condition for the velocity function in (4.2b), we find that $\tilde {F}_{1}$ should be varying linearly as

(4.11)\begin{equation} \tilde{F}_{1}(\tilde{\eta}_{1})={-}C_3 + C_4\tilde{\eta}_{1}. \end{equation}

This linear variation of $\tilde {F}_{1}$ for large $\tilde {\eta }_{1}$ can be seen in the inset of figure 8, which shows the variation of $\tilde {F}_{1}(\tilde {\eta }_{1})$ with $\tilde {\eta }_{1}$ obtained by solving (4.1). We obtain $C_3=1.04$ and $C_4=1.15$ from the $y$ intercept and the slope of the linear portion of this plot. Substituting (4.11) in (4.1c) and integrating gives us $\tilde {H}'_{1}(\tilde {\eta }_{1})$ at large $\tilde {\eta }_{1}$ as

(4.12)\begin{equation} \tilde{H}'(\tilde{\eta})=C_5\, {\rm e}^{{-}3C_4(\tilde{\eta}- C_3/C_4)^2/10}. \end{equation}

A curve fit to the values of $\tilde {H}'(\tilde {\eta })$ for $Pr \gg 1$, obtained by numerically solving (4.1), showed that $C_5=-0.38$. Converting (4.11) and (4.12) into $F(\eta )$ and $H'(\eta )$ by using (4.7) and (4.8), and then substituting in (3.11) and using (4.6), we get

(4.13)\begin{equation} \eta_{df}Pr^{1/5}=0.90+\sqrt{1.45W_0 ( 0.21(Th_mPr)^2 )}, \end{equation}

where $W_0(x)$ is the Lambert W function, which is the inverse function of $y\,{\rm e}^y=x$.

Equation (4.13) is the expression for $\eta _{df}$ in the large $\eta$ limit for $Pr \gg 1$. In figure 8, we compare (4.13) with $\eta _{df}$ obtained from (3.11) for $Pr=2$, $Pr=10$ and $Pr=100$, where the $F(\eta )$ and the $H'(\eta )$ corresponding to $Pr=2$, $Pr=10$ and $Pr=100$ used in (3.11) were obtained by numerically solving (2.28). As expected, the data of $Pr=100$ are very close to the curve of (4.13), which was derived for the limit of $Pr \gg 1$, while the data of $Pr=2$ deviate from (4.13). Thus, (4.13) seems to hold as an analytic expression for the dust-free region height for $Pr \ge 100$ for large $\eta _{df}$ values corresponding to higher $Th_m$. We now examine the scaling laws for the large and the small $\eta$ cases for low $Pr$.

4.2. Scaling for $Pr\ll 1$

4.2.1. Scaling for large $\eta$

For $Pr \ll 1$, Rotem & Claassen (Reference Rotem and Claassen1969) gives the outer boundary layer equations in terms of the dimensionless velocity $\tilde {\tilde {F_2}}$ in (D12a), pressure $\tilde {\tilde {G_2}}$ in (D13) and temperature function $\tilde {\tilde {H_2}}$ in (D12b) in the outer boundary layer as

(4.14a–c)\begin{equation} 3\tilde{\tilde{F_2}}\tilde{\tilde{F''_2}}-(\tilde{\tilde{F'_2}})^2= 2\tilde{\tilde{G_2}} - 2\tilde{\tilde{\eta_2}} \tilde{\tilde{G'_2}}, \quad \tilde{\tilde{H_2}} = \tilde{\tilde{G'_2}}, \quad \tilde{\tilde{H''_2}} + \frac{3}{5} \tilde{\tilde{F_2}}\tilde{\tilde{H'_2}}=0, \end{equation}

with the boundary conditions

(4.15a)\begin{equation} \tilde{\tilde{F_2}}=0, \quad \tilde{\tilde{H_2}}=\tilde{\tilde{G'_2}}=1 \quad \mathrm{at} \ \tilde{\tilde{\eta_2}} =0 \end{equation}

and

(4.15b)\begin{equation} \tilde{\tilde{F'_2}}=0, \quad \tilde{\tilde{H_2}}=\tilde{\tilde{G_2}}=0 \quad \mathrm{when} \ \tilde{\tilde{\eta_2}} \rightarrow \infty. \end{equation}

The subscript 2 in the above equations is for the outer boundary layer.

From (D15), (D3) and (D11c), we have

(4.16)\begin{equation} F'(\eta) = Pr^{{-}1/5}\tilde{\tilde{F'_2}}(\tilde{\tilde{\eta_2}}). \end{equation}

Similarly, we have

(4.17)\begin{equation} \eta=Pr^{{-}2/5}\tilde{\tilde{\eta_2}} \end{equation}

from (D11a), (D11b), (D12c) and (2.21). Using (4.17) in (4.16), integrating, and since $F(0)$ and $\tilde {\tilde {F_2}}(0)$ are equal to zero from (2.30a) and (4.15a), respectively, we get

(4.18)\begin{equation} F(\eta) = Pr^{{-}3/5}\tilde{\tilde{F_2}}(\tilde{\tilde{\eta_2}}). \end{equation}

Differentiating (D12b) and using (2.20) and (4.17), we get

(4.19)\begin{equation} H'(\eta)=Pr^{2/5}\tilde{\tilde{H'_2}}(\tilde{\tilde{\eta_2}}). \end{equation}

From the boundary condition for velocity function in (4.15b) for $\tilde {\tilde {\eta _2}} \rightarrow \infty$, we see that

(4.20)\begin{equation} \tilde{\tilde{F_2}} (\tilde{\tilde{\eta_2}})= C_6, \end{equation}

where the constant $C_6=1.78$ is obtained from the numerical solution of (4.14). Using (4.20) in (4.1c) and integrating gives

(4.21)\begin{equation} \tilde{\tilde{H'_2}}(\tilde{\tilde{\eta_2}}) = C_7\, {\rm e}^{{-}3 C_6\tilde{\tilde{\eta_2}}/5}, \end{equation}

where $C_7=-1.97$, as obtained by numerically solving (4.14) for $\tilde {\tilde {H'_2}}(\tilde {\tilde {\eta _2}})$.

Using (4.17), (4.18) and (4.19) to convert (4.20) and (4.21) to be in terms of $\eta$, and then substituting in (3.11), we get the dimensionless dust-free region height at large $\eta$ for $Pr \ll 1$ as

(4.22)\begin{equation} \eta_{df}=0.93Pr^{{-}2/5}\ln(1.84Th_mPr). \end{equation}

In figure 9(a), we compare (4.22) with the $\eta _{df}$ obtained from (3.11) by using the numerical solution of (2.28) for $Pr=0.1$ and $Pr=1$. Since (4.22) was obtained from the equations of Rotem & Claassen (Reference Rotem and Claassen1969) for $Pr \ll 1$ in (4.14), as expected, the numerical solution for $Pr=0.1$ agrees with the scaling relation (4.22) for large $\eta$. However, surprisingly, the data for $Pr=1$ also match with (4.22) for large $\eta$, suggesting that (4.22) is valid even for $Pr\le 1$. Thus, (4.22) gives the scaling for the dimensionless dust-free region height for $Pr\le 1$ for large $Th_m$.

Figure 9. Scaling of the dimensionless dust-free region height ($\eta _{df}$) at $Pr \ll 1$ in the large and the small $\eta$ limits. (a) Large $\eta$ limit. Numerical solutions at green $\circ$, $Pr=0.1$ and magenta $\diamond$, $Pr=1$; – – – –, scaling relation (4.22). (b) Small $\eta$ limit. Numerical solutions at green $\circ$, $Pr=0.1$ and magenta $\diamond$, $Pr=1$; – – – –, scaling relation (4.29). The inset shows the variation of $\tilde {\tilde{H}}^{\prime}_2(\tilde {\tilde {\eta }})$ with $\tilde {\tilde {\eta }}$ for $Pr\ll 1$, obtained by solving the outer boundary layer equations (4.14).

4.2.2. Scaling for small $\eta$

Rotem & Claassen (Reference Rotem and Claassen1969) provides the inner boundary layer equations for $Pr \ll 1$ as

(4.23a,b)\begin{equation} 5\tilde{\tilde{F^{\prime\prime\prime}_1}}+3\tilde{\tilde{F_1}}\tilde{\tilde{F^{\prime\prime}_1}}-\tilde{\tilde{F{^{\prime}}_1}}^2={-}1, \quad \tilde{\tilde{H_1}} = 1, \end{equation}

subjected to the boundary conditions,

(4.24a)\begin{equation} \tilde{\tilde{F_1}}=\tilde{\tilde{F'_1}}=0 \quad \mathrm{at} \ \tilde{\tilde{\eta_1}} =0 \end{equation}

and

(4.24b)\begin{equation} \tilde{\tilde{F'_1}}=1 \quad \mathrm{when} \ \tilde{\tilde{\eta_1}} \rightarrow \infty. \end{equation}

Subscript $1$ in the above equations denote the inner boundary layer. However, (4.23b) implies that the temperature gradient inside the inner velocity boundary layer, $\tilde {\tilde {H'_1}}=0$; this is quite unphysical since the maximum temperature gradient in natural convection flows occurs at the wall. So we assume that the temperature gradient in the outer velocity boundary layer, which approaches a constant value at lower $\tilde {\tilde {\eta _2}}$, is maintained throughout the inner thermal boundary layer. By numerically solving the set of outer boundary layer equations (4.14), we get

(4.25)\begin{equation} \tilde{\tilde{H'_2}}(\tilde{\tilde{\eta_2}})=C_8={-}0.57, \end{equation}

as $\tilde {\tilde {\eta _2}} \rightarrow 0$, shown as the $y$-intercept in the inset of figure 9(b). From a local analysis of (4.23) near $\tilde {\tilde {\eta _1}}=0$, given in Appendix C.3, we see that $\tilde {\tilde {F_1}}$ can be approximated to

(4.26)\begin{equation} \tilde{\tilde{F_1}}(\tilde{\tilde{\eta_1}})=C_9\tilde{\tilde{\eta_1}} ^{2}. \end{equation}

As given by Rotem & Claassen (Reference Rotem and Claassen1969),

(4.27)\begin{equation} F(\eta) = Pr^{{-}1/10}\sqrt{\tilde{\tilde{F'_2}}(0)}\,\tilde{\tilde{F_1}}(\tilde{\tilde{\eta_1}}) \end{equation}

and

(4.28)\begin{equation} \eta=\frac{Pr^{1/10}\tilde{\tilde{\eta_1}}}{\sqrt{\tilde{\tilde{F'_2}}(0)}}. \end{equation}

Using (4.27) to convert (4.26) into $F(\eta )$, (4.19) to convert (4.25) into $H'(\eta )$, substituting these in (3.11), and using (4.28), we get

(4.29)\begin{equation} \eta_{df}=1.75Pr^{7/20}Th_m^{1/2} \end{equation}

as the scaling law to predict $\eta _{df}$ for $Pr \ll 1$ in the small $\eta$ limit. Here, the pre-factor $1.75$ is chosen to match (4.29) with the data from the solution of (3.11) for $Pr<1$, plotted in figure 9(b). Figure 9(b) shows the comparison of (4.29) with the dust-free region height obtained from the numerical solutions of (3.11) for $Pr=0.1$ and $Pr=1$. We again note from the figure that (4.29) is valid for any $Pr\le 1$, even though it was derived for $Pr \ll 1$.

4.2.3. Large $\eta$ scaling for intermediate $Pr$ ($1< Pr<10$)

We saw in §§ 4.1.1 and 4.2.2 that the scaling laws for the dust-free region heights obtained in the small $\eta$ limit for $Pr \gg 1$ in (4.9) and for $Pr \ll 1$ in (4.29) are valid for any $Pr>1$ and $Pr<1$. Similarly, the scaling law for the large $\eta$ limit for $Pr \ll 1$ in (4.22) was shown to be valid for any $Pr<1$. However, the scaling law in the large $\eta$ limit for $Pr \gg 1$ in (4.13) deviates from the simulation data for $Pr=2$ (see figure 8). Then in the intermediate range of $Pr$ values which are not much larger than one, for large $Th_m$ that push the $\eta _{df}$ to large $\eta$ limits, the prediction of $\eta _{df}$ becomes difficult with the results presented so far. To overcome this limitation, we now obtain an expression, by trial and error using numerical solutions of $F(\eta )$ and $H'(\eta )$ obtained from (2.28), which collapse the $\eta _{df}$ data from the solutions of (3.11) for $1< Pr<10$. A logarithmic function,

(4.30)\begin{equation} \eta_{df} Pr^{0.39} = 0.68 \ln(41.2 Th_m Pr^{1.5}), \end{equation}

and the corresponding data from the solution of (3.11) are compared in figure 10 for the large $\eta$ limit for $2< Pr<10$; the data of simulations collapse on to (4.30) fairly well for this range of intermediate $Pr$.

Figure 10. Comparison of the scaling of the dimensionless dust-free region height for $1< Pr<10$ with $\eta _{df}$ obtained from the numerical solution; blue $\triangle$, $Pr=2$; red $\circ$, $Pr=5$; green $*$, $Pr=7$; magenta $\diamond$, $Pr=10$; – – – –, (4.30).