3.1. Time-averaged flow structures

A study by Hartlep et al. (Reference Hartlep, Tilgner and Busse2005) reveals the spatial and temporal changes of structures in RBC for Rayleigh numbers up to $10^{7}$ and Prandtl number in the range of $0.7$ and $60$ for large-aspect-ratio cells. In particular, when $Pr\approx 1$ , straight rolls occupy the domain at the onset of convection. Subsequently, the regularly organized roll-like structures are replaced by a random array of transient plumes as the Rayleigh number increases. When the Rayleigh number increases beyond $10^{4}$ , the flow becomes chaotic. In our study, at parameters $Ra=10^{6}$ , $Pr=0.7$ and $\varGamma =6$ , the flow regime was soft turbulence, in which the flow was chaotic and the thermal plumes moved randomly. As mentioned by Hartlep, Tilgner & Busse (Reference Hartlep, Tilgner and Busse2003), a long-term average of the fluid velocities converges to zero because the flow runs in different directions at different times. Typical instantaneous temperature distributions caused by thermal plumes are shown in figures 1(a) and 1(b) in terms of isosurfaces and two-dimensional distribution at the central plane, respectively. The red structures drawn on the isosurface at $T/\Delta T=0.8$ indicate that hot plumes converge horizontally in a sheet-like pattern near the bottom plate and create a fine network across the plate, ejecting thermal plumes vertically, as explained by Zocchi et al. (Reference Zocchi, Moses and Libchaber1990). Then detached thermal plumes accelerated vertically and crossed the bulk region in the shape of a mushroom. Finally, they impinge on the top plate and spread out horizontally. Owing to the symmetric nature of RBC, cold plumes behave similarly in opposite directions. Therefore, the distribution of hot rising and cold falling fluids is symmetric in the vertical direction, as shown in figures 1(a) and 1(b), and this symmetry is referred to as the up–down symmetry.

However, when particles are laden in the flow, we find that a polygonal cellular pattern emerges, as shown in the examples given in figures 1(c) and 1(d). The particle size and mass loading were $d_{p}=0.003D$ and $\varPhi _{m}=0.12$ . In this example, the observed pattern consists of square cells in which hot fluid soars up in the centre of each cell and cold fluid sinks in the peripheral region of the cell, breaking the up–down symmetry. This cellular pattern persisted over time. This stationary pattern was caused by the settling of the particles. The occurrence of a hexagonal pattern under symmetric boundary conditions was experimentally observed by Berdnikov & Kirdiashkin (Reference Berdnikov and Kirdiashkin1979). In their experiment, aluminium particles of the size of 10 microns were added for visualization to the 4 mm ethyl alcohol layer with an aspect ratio larger than $30$ at Rayleigh number of $10^{4}$ . They observed that the agitated solid particles, while settling, altered the flow pattern to a hexagonal cell structure, with a descending flow in the centre of each cell. Although their Rayleigh number is much lower than ours of $10^6$ and their Prandtl number is much higher than ours, their hexagonal cells are the only examples of polygonal patterns caused by settling particles reported thus far.

Figure 1. Instantaneous temperature distribution: (a,b) of a particle-free flow; (c,d) of a particle-laden flow at $Ra=10^6$ when $\varPhi _{m}=0.12$ , $d_{p}/D=0.003$ . Isosurface at $T/\Delta T=0.8$ with cross-sectional distribution at selected vertical planes are shown in (a) and (c). Temperature distribution at the central plane of $y/D=0.5$ is shown in (b) and (d).

Figure 2. Isosurfaces of time-averaged temperature at $Ra=10^6$ . Hot thermal plumes are drawn in red colour at $T/\Delta T = 0.8$ and cold ones are in blue colour at $T/\Delta T = 0.2$ . The two isosurfaces are more separately drawn than real distance for better visualization.

Figure 3. Contours of time-averaged vertical velocity at an $xz$ -plane of $y/D=0.5$ at $Ra=10^6$ . The red area is where the rising flow occurs, and the blue area is where the falling flow occurs. The range of contours is between $-0.2U_{\!f}$ (blue) and $0.2U_{\!f}$ (red).

The polygonal cellular patterns observed in our simulations for the three particle sizes and three mass loadings at $Ra=10^6$ are illustrated on the isosurfaces of the time-averaged temperature and two-dimensional distribution of the time-averaged vertical velocity at the half-channel height, as shown in figures 2 and 3, respectively. Time averaging for approximately $1200D/U_{\!f}$ yielded a smooth distribution of isosurfaces at $T/\Delta T=0.2$ (blue) and 0.8 (red). When a polygonal cell pattern was observed for $d_p/D = 0.0015$ (small) and 0.003 (intermediate), the hot plumes were always concentrated at the centre of each cell, whereas the cold plumes were distributed on the periphery of each cell, as shown in figure 2. When roll structures were observed for $d_p/D= 0.0045$ (large), straight rolls were observed for $\varPhi _m=0.06$ whereas crooked rolls and straight but knotted rolls were obtained for $\varPhi _m=0.12$ and 0.18, respectively. From the vertical velocity distribution at the channel centre height in figure 3, the shape of the cells could be better identified, and a polygonal structure was observed for cases with small and intermediate particles ( $d_p/D=0.0015$ and 0.003), whereas a roll structure was found for large particles ( $d_p/D=0.0045$ ). For all particle sizes, the sizes of the cells or rolls tended to decrease as the mass loading increased. The polygonal cell structures indicated that the up–down symmetry was broken, whereas in the roll structures, the symmetry was maintained. For $Ra=10^7$ and $10^8$ , similar polygonal cell structures were observed as illustrated in Appendix A, confirming that this kind of polygonal cell structures is not a phenomenon observable only for low Rayleigh numbers.

To quantify the size of the observed cells or rolls, we introduce a length scale $\lambda$ and a corresponding wavenumber $k$ for the cell structures, defined as

(3.1) \begin{equation} \lambda = \left ( \frac {L_x L_z}{2n} \right )^{1/2},\quad k =\frac {2\pi }{\lambda }, \end{equation}

where $n$ is the number of polygonal cells in the entire domain. For roll structures, $\lambda$ is defined as the horizontal wavelength, which is the lateral size of a pair of rolls. The values of $n, \lambda$ and $k$ for all cases considered in this study are listed in table 3. For reference, the corresponding values of $\lambda$ and $k$ for the most unstable mode of unladen RBC at the critical Rayleigh number $Ra_{c} (=1708)$ from linear stability analysis are $\lambda _c/D=2.016$ and $k_c D=3.117$ , respectively (Drazin & Reid Reference Drazin and Reid1981). To get some clues on why the settling particles induce such cell or roll structures, we carried out the linear stability analysis of a particle-laden RBC as provided in Appendix B. Because the particle Stokes number $St (=\tau _p \nu /D^2)$ is much smaller than one as shown in table 3, the so-called flow of particles exists, guaranteeing a smooth distribution of the particle velocity (Fouxon et al. Reference Fouxon, Park, Harduf and Lee2015). This allows the particle feedback force that is the Stokes drag force to be explicitly expressed in terms of the local slip velocity field, as discussed in Appendix B. We found that the main parameter quantifying the effect of particles on fluid motion in the linear stability analysis was $\phi _m (=\varPhi _m S)$ , where $S$ denotes the non-dimensionalized particle settling velocity $\tau ^{\prime}_p g D/\nu$ , which is listed in table 3. The stability analysis indicates that the critical Rayleigh number increases with $\phi _m$ as shown in Appendix B. Therefore, settling particles tend to stabilize the flow, which is consistent with the previous findings by Prakhar & Prosperetti (Reference Prakhar and Prosperetti2021) and Raza et al. (Reference Raza, Hirata and Calzavarini2024). For the range $\phi _m = 200 - 2000$ considered in the study, $Ra_c \simeq 35\,000 - 350\,000$ . The Rayleigh number of our cases, however, is $Ra=10^6$ , which is not much higher than $Ra_c$ , suggesting that laminar cells or roll structures are preferably formed because the turbulence is expected to be weak. However, the critical wavenumber predicted by the stability analysis is relatively constant at $k_c D \simeq 3.7$ or $\lambda _c/D\simeq 1.7$ for the considered range of $\phi _m$ as shown in table 7 and figure 18(b), which cannot capture the behaviour of $k D$ varying with the mass loading listed in table 3. Another limitation of the linear stability analysis is the indeterminacy of the cell shape. Despite these limitations, the linear stability analysis suggests that the settling particles stabilize the flow and thus suppress turbulence, possibly leading to the formation of cell or roll structures. The value of $\lambda _c/D (=1.7)$ obtained from the linear stability analysis is also compared in the figure 4, indicating that a comparable $\lambda /D$ is observed only for $d_p/D=0.0015$ and 0.003 at $\varPhi _m=0.18$ for which the total number of particles is large. It can be conjectured that the discrete nature of the particle feedback force, which cannot be reflected in the linear stability analysis, is responsible for the larger cell or roll sizes compared with those predicted by the linear stability analysis.

Figure 4. The cell size $\lambda /D$ shown in table 3 as a function of $\varPhi _m$ at $Ra=10^6$ . Fitting lines, $\lambda /D \sim 0.5 \varPhi _m^{-2/3}$ (solid line) for the cell structure and $\lambda /D \sim 1.5 \varPhi _m^{-1/3}$ (dashed line) for the roll structure are provided to guide the decreasing trend. The value of $\lambda /D =1.7$ proposed by the linear stability analysis is shown together in the dash–dot line.

Figure 5. Statistics of time-averaged field when mass loading $\varPhi _{m}=0.06$ , $0.12$ and $0.18$ at three different particle sizes at $Ra=10^6$ . (a) Volume-averaged temperature in the whole flow domain, (b) horizontally averaged temperature, (c) horizontally averaged horizontal components of velocity, (d) horizontally averaged vertical component of velocity. Red, blue and green coloured symbols or lines denote data for small ( $d_p/D=0.0015$ ), intermediate ( $d_p/D=0.003$ ) and large ( $d_p/D=0.0045$ ) sized particles, respectively. Solid, dashed and dotted lines for each colour indicate mass loadings of 0.06, 0.12 and 0.18, respectively. These legends apply to all figures throughout the paper.

For the quantitative analysis of the shape of cells or rolls, figures 5(a) and 5(b) provide the statistics of the time-averaged temperature $\overline {T}$ defined by

(3.2) \begin{equation} \overline {T}(x,y,z) = \frac {1}{\mathcal {T}} \int _0^{\mathcal {T}} T(x,y,z,t)\, {\rm d}t \end{equation}

from which the horizontal and volume averages are defined by

(3.3) \begin{equation} \langle \overline {T}\rangle _{A}(y) = \frac {1}{L_x L_z} \int _0^{L_z} \int _0^{L_x} \overline {T}(x,y,z)\,{\rm d}x\,{\rm d}z ,\,\, \langle \overline {T}\rangle _{V} = \frac {1}{D} \int _0^{D} \langle \overline {T}\rangle _A(y)\,{\rm d}y . \end{equation}

Throughout this paper, the overline and bracket represent the time-and space-averaged operations, respectively. The values of the volume-averaged and horizontally averaged temperatures are shown in figures 5(a) and 5(b) clearly indicate a bias of $\overline {T}$ towards the hot bottom temperature. The volume-averaged temperature for all cases exhibiting cell or roll structures is greater than 0.5 and the horizontal average is greater than 0.5 in most regions except for a small region near the top wall, breaking the up–down symmetry. The asymmetry is more pronounced for cases showing the polygonal cell structure for $d_p/D=0.0015$ and 0.003 because hot plumes are concentrated in the centre of each cell and cold plumes are distributed along the edge of each cell, as shown in figure 2. This also confirms that the hot plumes are stronger than the cold plumes, resulting in a volume-averaged temperature greater than 0.5. Furthermore, smaller particles induced stronger hot plumes. As more cells were formed with an increase in mass loading, the asymmetry became stronger. However, when the roll structure was observed for $d_p/D=0.0045$ , a relatively symmetric temperature distribution was observed.

Figures 5(c) and 5(d) show the distribution of the horizontally averaged statistics of the time-averaged velocity field for the cell and roll structures. Because no mean flow was caused by the particles ( $\langle \overline {u_i}\rangle _{V} = \langle \overline {u_i}\rangle _A=0$ ), the standard deviation of the velocity components was investigated. The distributions of the standard deviations of the horizontal and vertical velocity components shown in figures 5(c) and 5(d) clearly indicate active horizontal motion near the walls and mainly vertical motion in the core, which is a typical behaviour of the cell or roll structure. A symmetric distribution was maintained for the roll structure ( $d_p/D=0.0045$ ), whereas an asymmetric distribution was observed for the cell structures ( $d_p/D=0.0015$ and 0.003). For $d_p/D=0.0015$ , the positions exhibiting the maximum vertical velocity and the minimum horizontal velocity were found slightly below the channel centre. Although the cases for $d_p/D=0.0015$ and 0.003 show similar cell structures at the same mass loading, as shown in figure 3, the statistics of the standard deviation show different behaviours, which might be due to detailed differences in the cell structure, such as the shape of the cross-section of hot plumes. For $d_p/D=0.0015$ , hot plumes in the centre of a cell have a circular cross-section, but for $d_p/D=0.003$ , the cross-section of the hot plume has a cross shape, which is similar to the shape of a cold plume, as shown in figure 3.

3.2. Turbulence

In this section, we investigate the turbulence behaviour for cell or roll structures defined by the fluctuations in the time-averaged fields through decomposition,

(3.4) \begin{equation} u_i(x,y,z,t) = \overline {u}_i(x,y,z) + u_i^{\prime}(x,y,z,t), \quad T(x,y,z,t) = \overline {T}(x,y,z)+T'(x,y,z,t), \end{equation}

where $\overline {u}_i(x,y,z)$ and $\overline {T}(x,y,z)$ correspond to the fields of the stationary cells or roll structures. In the particle-free case, $\overline {u}_i=0$ . The horizontal and temporally averaged turbulence intensities for the cell and roll structures are shown in figure 6 with those for the particle-free case. For the particle-free case, the intensities in the horizontal and vertical directions were comparable because the randomly occurring hot or cold plumes and distributions were such that the horizontal components peaked near the top and bottom walls. Conversely, the vertical component exhibited a maximum value in the channel centre. The intensities of the cell or roll structures demonstrate that the horizontal components of the velocity are significantly suppressed by the settling particles, whereas the vertical component of the velocity is moderately suppressed or slightly enhanced, depending on the particle size. This observation indicates that hot or cold plumes randomly occurring in the particle-free case are suppressed by the settling particles. Instead, steady plume activity is caused by particles in the form of a cell or roll structure. The relatively less suppressed vertical component is probably due to the discrete distribution of the particle feedback forces. The nearly uniform distribution of the vertical component across the channel and the monotonic increase in the particle size support this conjecture. Furthermore, it is noticeable that all three intensity components are stronger near the bottom wall than near the top wall for the cell structures for $d_p/D=0.0015$ and 0.003, whereas they are approximately symmetric for the roll structure for $d_p/D=0.0045$ . By contrast, the temperature fluctuation is not substantially suppressed by the particles, as shown in figure 6(c). The relatively slight changes in mean temperature shown in figure 5 are responsible for this behaviour. The temperature fluctuation near the bottom wall for $d_p/D=0.0015$ is suppressed more than that near the top wall, which is the opposite trend to the behaviour of the turbulence intensities shown in figures 6(a) and 6(b). A strong asymmetry is observed in the behaviour of the temperature skewness for $d_p/D=0.0015$ as shown in figure 6(d). Overall, for the roll structure with $d_p/D=0.0045$ , most turbulence statistics remain symmetric, whereas for the cell structure with $d_p/D=0.0015$ and 0.003, the turbulence is not symmetric. The observation of turbulence in the cell structure combined with the time-averaged fields indicated that both the time-averaged flow and turbulence near the bottom wall were more intense than those near the top wall, whereas both the time-averaged temperature and fluctuation near the bottom wall were less intense than those near the top wall.

Figure 6. Standard deviations of fluctuating velocity and temperature when mass loading $\varPhi _{m}=0.06$ , $0.12$ and $0.18$ at different particle sizes at $Ra=10^6$ : (a,b) for velocity components; (c) for temperature; (d) skewness factor of temperature, respectively.

Figure 7. One-dimensional spectra of fluctuating velocity and temperature at the middle plane $y/D=0.5$ when mass loading $\varPhi _{m}=0.06$ , $0.12$ and $0.18$ at different particle sizes at $Ra=10^6$ : (a,b,c) for velocity components; (d) for temperature, respectively.

The scale characteristics of the turbulence were investigated through one-dimensional spectra of the turbulent velocity components and temperature fluctuations on the horizontal plane at $y/D=0.5$ as shown in figure 7. One-dimensional spectra were defined as twice the one-dimensional Fourier transform of the two-point correlation,

(3.5) \begin{equation} E_{ij}(\kappa _{x}) = \frac {1}{\pi } \int _{-\infty }^{\infty } R_{ij}(r_{x})e^{-i\kappa _{x}r_{x}}{\rm d}r_{x},\quad E_{T'}(\kappa _x) = \frac {1}{\pi } \int _{-\infty }^{\infty } R_{T'}(r_{x})e^{-i\kappa _{x}r_{x}}{\rm d}r_{x}, \end{equation}

where $R_{ij}$ and $R_{T'}$ are two-point correlations,

(3.6) \begin{equation} R_{ij}(r_{x}) = \langle u_{i}^{\prime}(\boldsymbol{x},t)u'_{\!j}(\boldsymbol{x}+\boldsymbol{e}_{x}r_{x},t) \rangle _{A,t},\,\, R_{T'}(r_{x}) = \langle T'(\boldsymbol{x},t)T'(\boldsymbol{x}+\boldsymbol{e}_{x}r_{x},t) \rangle _{A,t}, \end{equation}

where $\kappa _x$ denotes the wavenumber in $x\hbox{-}$ direction. $E_{11}=E_{u'}, E_{22}=E_{v'}$ and $E_{33}=E_{w'}$ in figure 7. From a comparison of the energy spectra with those for the particle-free case, it is easily recognized that the settling particles suppress large-scale motions while enhancing small-scale motions in all cases. In particular, small-scale vertical motion is substantially increased owing to the quick settling of heavy particles. This trend became more pronounced as the particle size and mass loading increased. However, the large-scale horizontal motion tended to be more suppressed with mass loading. We conjecture that the modification of the vertical motion by the particles is rather direct and local because of the small size of the particles, whereas the horizontal motion is affected by the modified pressure to satisfy the divergence-free constraint, which is global. However, the spectrum of the temperature fluctuations was not modified significantly by the settling particles, as shown in figure 7(c). Large-scale fluctuations were slightly suppressed in all the cases.

3.3. Heat flux

In this section, we investigate the behaviour of the heat flux modified by the settling particles through mechanical coupling only, because we did not consider thermal coupling in this study. From (2.3), the heat flux vector $H_j$ is defined as follows:

(3.7) \begin{equation} \frac {\partial T}{\partial t} + \frac {\partial H_j}{\partial x_{\!j}} = 0, \quad H_j = u_{\!j} T -\kappa \frac {\partial T}{\partial x_{\!j}}. \end{equation}

In particular, the vertical heat flux is important because it determines the Nusselt number,

(3.8) \begin{equation} Nu = \frac { \langle \overline {H}_2 \rangle _{A} }{\kappa \Delta T/D} = \sqrt {Ra Pr} \frac { \langle \overline {H}_2 \rangle _{A} }{U_f \Delta T} , \end{equation}

where

(3.9) \begin{equation} \langle \overline {H}_2 \rangle _{A} (y) = \langle \overline {vT} \rangle _A (y) - \kappa \frac {{\rm d} \langle \overline {T} \rangle _{A}}{{\rm d}y} (y) = - \kappa \left . \frac {{\rm d}\langle \overline {T} \rangle _{A}}{{\rm d}y} \right |_{y=0} = - \kappa \left . \frac {{\rm d}\langle \overline {T} \rangle _{A}}{{\rm d}y} \right |_{y=D}, \end{equation}

where the last two equalities hold under stationary conditions. The behaviour of $Nu$ modified by the settling particles is shown in figure 8 and listed in table 3. The settling particles affected the heat flux and augmented $Nu$ for $d_p/D=0.003$ and 0.0045, whereas they suppressed $Nu$ for $d_p/D=0.0015$ compared with that for the particle-free case. Given that thermal coupling was not considered, the cell or roll structures formed by the particles modified $Nu$ significantly through momentum coupling. Roll structures seem to enhance heat transfer more effectively than cell structures. However, similar cell structures affect heat transfer differently depending on the size of the particles: intermediate-sized particles ( $d_p/D=0.003$ ) enhance heat transfer with an increase in mass loading, whereas small particles ( $d_p/D=0.0015$ ) suppress heat transfer. A similar decrease in $Nu$ by laden small particles was reported in Oresta & Prosperetti (Reference Oresta and Prosperetti2013) and Park et al. (Reference Park, O’Keefe and Richter2018), but the basic mechanism behind the similar phenomenon is different because no cell or roll structure was observed in their studies.

Figure 8. Nusselt number distribution with mass loading for different particle diameters at $Ra=10^6$ . The letter near each symbol denotes cell structure: ‘R’ for roll structure, ‘S’ for square cell, ‘P’ for pentagonal cell and ‘H’ for hexagonal cell structures. Nusselt number for the case without particles, 8.14 is shown for comparison.

Figure 9. Time-averaged local heat flux at boundary walls for different particle sizes with fixed mass loading of $\varPhi _{m}=0.06$ , (a,b) for $d_{p}/D=0.0015$ , (c,d) for $d_{p}/D=0.003$ , (e,f) for $d_{p}/D=0.0045$ at $Ra=10^6$ ; (a,c,e) bottom wall, (b,d,f) top wall.

Figure 10. Horizontally averaged local convective heat flux at $Ra=10^6$ : (a) total convective heat flux, (b) conductive heat flux, (c) heat flux of the time-averaged flow, (d) turbulent heat flux.

The contributions of the cell or roll structures formed by the settling particles to the vertical heat transfer can be investigated in terms of the time-averaged heat flux distribution at the walls, as shown in figure 9. The edges of the hexagonal or square cells were identified by the region showing a locally high heat flux at the bottom wall (figures 9 a and 9 c) and locally low heat flux at the top wall (figures 9 b and 9 d). At the centre of the cells, a locally low heat flux at the bottom wall and a locally high heat flux at the top wall were observed. This suggests that hot plumes contribute strongly to the heat transfer at the top wall in the centre of each cell, whereas cold plumes enhance the heat transfer at the bottom wall along the edge of each cell. However, when the roll structures are formed, such asymmetry between the hot and cold plumes is not observed in the cell structures, as shown in figures 9(e) and 9( f). Furthermore, it appears that the roll structures are more efficient at transporting heat than the cell structures under the same mass loading of particles.

For a more detailed investigation of the heat transport behaviour, the contribution of the turbulent heat flux was separated from the total heat flux through the following decomposition:

(3.10) \begin{equation} \langle \overline {H}_2 \rangle _A = \langle \overline {vT} \rangle _A - \kappa \frac {{\rm d}\langle \overline {T} \rangle _A }{{\rm d}y} = \langle \overline {v}\overline {T} \rangle _A + \langle \overline {v'T'} \rangle _A - \kappa \frac {{\rm d}\langle \overline {T} \rangle _A }{{\rm d}y} . \end{equation}

The distributions of the normalized components of heat flux are shown in figure 10. The convective transport of heat in the core region (figure 10 a) determines the heat flux at both boundaries or, equivalently, the total heat flux (figure 10 b). Notably, the convective heat transport increases with mass loading for particles with $d_p/D=0.003$ and 0.0045, whereas it decreases slightly with mass loading for particles with $d_p/D=0.0015$ . More detailed behaviour can be investigated in the decomposed elements of the convective heat transport, the heat transport due to the cell or roll structures $\langle \overline {v} \overline {T} \rangle _A$ and the turbulent heat transport $\langle \overline {v'T'} \rangle _A$ shown in figures 10(c) and 10(d), respectively. The heat transport due to the cell or roll structures was enhanced with mass loading for all particle sizes, whereas turbulent heat transport was not augmented with mass loading. In the case of particles with $d_p/D=0.0015$ , the turbulent transport was suppressed even with an increase in mass loading, which was responsible for the suppression of the total heat flux. This is obviously due to the suppression of the turbulence of both the vertical momentum and temperature with an increase in the mass loading for $d_p/D=0.0015$ shown in figures 6(b) and 6(c), unlike the cases for $d_p/D=0.003$ and 0.0045. Therefore, the decrease in Nusselt number for $d_p/D=0.0015$ with the mass loading reported in figure 8 is mainly attributed to the suppression of turbulence by the smallest particles, compensating for the augmented transport by the cell or roll structures.

3.4. Formation of polygonal structures

In this section, we investigate the formation of polygonal cell structures resulting from the particle feedbacks. We focused on the initial temporal behaviour of the emerging flow structures owing to the activation of feedback forces induced by sedimenting particles. Figure 11 provides an example of the sequential representation of the flow structures in the midplane temperature distribution, beginning at $tD/u_{\!f}=0$ , coinciding with the onset of the feedback force action and extending through $tD/u_{\!f}=300$ , at which point the manifestation of a square-cell structure becomes apparent. It was noticed that randomly distributed hot and cold plumes in the particle-free flow were disturbed by settling particles in the early period, but they became stronger. Soon after, the hot plumes merged, while the cold plumes aligned and organized into a square cell structure. The formed square cells resembled the initial plume distribution, suggesting that the initial field was responsible for determining the phases of the formed square cells. Another test with a different initial field yielded similar square cell structures with different phases.

Figure 11. Temperature distribution at midplane at different times during formation of cell structures when $d_{p}/D=0.003$ and $\varPhi _{m}=0.06$ at $Ra=10^6$ : (a) $tD/U_{\!f}=0$ , (b) $tD/U_{\!f}=4$ , (c) $tD/U_{\!f}=10$ , (d) $tD/U_{\!f}=39$ , (e) $tD/U_{\!f}=100$ , ( f) $tD/U_{\!f}=300$ .

Figure 12. Transient behaviour of (a) magnitude of the horizontal components of velocity $(u^2 + w^2)/2$ , (b) magnitude of the wall-normal velocity $v^2$ , (c) Nusselt number, (d) volume-averaged enstrophy at $Ra=10^6$ . The constant lines at the later time of each figure are the flow component calculated from the time-averaged field in the steady state.

The transient behaviours of the mean kinetic energy, Nusselt number and enstrophy are demonstrated in figure 12 to provide an idea of how particles interact with the fluid, leading to the formation of cell or roll structures. The volume-averaged horizontal kinetic energy of the particle-free initial flow is quickly suppressed by the settling particles within $tD/U_f=5$ and then slowly recovers slightly, retaining lower levels than those of the particle-free flow for all cases, as shown in figure 12(a). The flat lines in the lower and right-hand corners indicate the mean kinetic energies of the time-averaged flows. A comparison of the levels of the total kinetic energy and the kinetic energy of the time-averaged flow indicates that the turbulent kinetic energy is relatively small and that the level of turbulent kinetic energy increases with mass loading for all particle sizes. However, the behaviour of the mean vertical kinetic energy was significantly different from that of the horizontal kinetic energy, as shown in figure 12(b). For $d_p/D=0.0015$ and 0.003, the vertical kinetic energy was initially suppressed substantially and slightly, respectively, whereas for $d_p/D=0.0045$ , the vertical kinetic energy increased. For all particle sizes, the vertical kinetic energy peaked at approximately $tD/U_f=10$ and then slowly decreased, converging to a stationary value. From the comparison with the mean vertical kinetic energy of the time-averaged flow, the turbulent kinetic energy of the converged flow is non-negligible owing to the direct interaction between the settling particles and the fluid. Compared with the vertical kinetic energy, however, the horizontal kinetic energy is substantially suppressed by the settling particles, making the flow more laminar. Then, the resulting flow becomes unstable, giving way to cellular convection as suggested by the linear stability analysis. However, the linear stability analysis does not determine the type of cell structure.

The Nusselt number initially dropped quickly to the minimum value at $tD/U_f=3$ and then suddenly increased to the maximum value at $tD/U_f=10$ , eventually converging to a stationary value for all particle sizes, as shown in figure 12(c). The peak values of the Nusselt number and vertical kinetic energy are observed at the same instance of $tD/U_f=10$ , confirming that the vertical motion of the fluid is critical for the vertical transport of heat. The mean enstrophy was significantly increased by particles with sizes of $d_p/D=0.003$ and 0.0045, whereas it was slightly decreased by particles with $d_p/D=0.0015$ as shown in figure 12(d).

Figure 13. Thermal plume structures with effective feedback force vectors in the cross-section of plumes at $Ra=10^6$ : (a) hot plume surface at $T/\Delta T=0.85$ and a cross-section and (b) cold plume surface at $T/\Delta T=0.15$ and a cross-section, for $d_{p}/D=0.003$ and $\varPhi _{m}=0.06$ (square cell structures); (c) hot plume surface and a cross-section and (d) cold plume surface and a cross-section for $d_p/D=0.0045$ and $\varPhi _m=0.06$ (roll structures). A thick solid line in (a) and (c) denotes a contour of $T/\Delta T=0.85$ and that in (b) and (d) a contour of $T/\Delta T=0.15$ .

To understand how the cell structure was formed by the settling particles, the effect of the feedback force was investigated. The vertical component of the feedback force plays a critical role in the flow modification. The effective feedback force in the vertical direction defined by $f^{\prime\prime}_2( \equiv f_2 - \langle f_2 \rangle _A )$ directly affects the flow because the horizontally averaged feedback force $\langle f_2 \rangle _A$ owing to the fast-settling particles, contributes to the buildup of the vertical mean pressure gradient owing to the confined geometry. An example distribution of the vertical effective feedback force in terms of the negative drag force acting on particles $-(F^k_{D,2} - \rho \Delta V \langle f_2 \rangle _A /m_p )$ in the cross-sections of hot and cold plumes for the case showing square cell structure ( $d_p/D=0.003, \varPhi _m=0.06$ ) is demonstrated in figures 13(a) and 13(b), where $\Delta V$ is the volume of the numerical cell. The first noticeable difference in the effective feedback force between the hot and cold plumes is in the opposite direction of the forces, which is caused by the fluid motion associated with the hot and cold plumes. When the particles are released at the top wall, where no vertical motion of the fluid is observed, the initial settling velocity of the particles is nearly the same as the terminal velocity in the still fluid. However, particles soon undergo rising and descending fluid motions in the hot and cold plume regions, respectively. Therefore, particles in the hot plume experience more drag because of the increased velocity difference, whereas particles in the cold plume experience less drag because of the decreased velocity difference, resulting in a downward effective feedback force in the hot plume and an upward effective feedback force in the cold plume in most regions, as shown in figures 13(a) and 13(b). However, this situation is reversed when the particles reach the near-bottom region, where the vertical fluid motion is suddenly suppressed by the bottom wall. Therefore, the vertical effective feedback force in the hot plume near the bottom wall is upward, whereas that in the cold plume near the wall is downward, as shown in figures 13(a) and 13(b). Interestingly, the upward force in the hot plume near the bottom wall enhanced the hot plume activity, whereas the downward force in the cold plume near the bottom wall did not affect any activity. This difference might help explain the more concentrated and stronger vertical motion of the hot plume and the widely distributed and weak vertical motion of the cold plume typically found in cell structures. For the roll structures, it is observed that such a difference in the feedback force near the bottom wall is not strong enough to break the symmetry, because the roll structures are found in cases with larger particles with large inertia settling so quickly that they are hardly observed near the bottom wall, as shown in figures 13(c) and 13(d). This suggests that small particles, which have small inertia and thus are more likely affected by the buoyancy force, tend to form polygonal cell structures through the asymmetric feedback forces near the wall, while large particles overcome the buoyancy force due strong inertia and create the symmetric feedback force, yielding symmetric structures such as the roll structure.

3.5. Energy budget analysis

Finally, to further analyse the interaction between the particles and fluid, we investigated the energy exchange between the settling particles and fluid through an energy budget analysis. As we did not consider the thermal interactions between them, we focused on the kinetic energy budget. The kinetic energy of the fluid can be obtained by multiplying $u_i$ by (2.2) and using (2.1),

(3.11) \begin{equation} \frac {\partial } {\partial t} \!\left (\frac {u_i u_i}{2}\!\right )+\frac {\partial }{\partial x_{\!j}}\!\left ( u_{\!j} \frac {u_i u_i}{2}\!\right ) = -\frac {1}{\rho } \frac {\partial (u_i p)}{\partial x_{i}}+\nu \frac {\partial }{\partial x_{\!j}} \left (\! u_i \frac {\partial u_{i}}{\partial x_{\!j}} \!\right ) - \nu \frac {\partial u_i}{\partial x_{\!j}} \frac {\partial u_i}{\partial x_{\!j}}+g \beta u_2 T +u_i f_{i}, \end{equation}

which can be averaged over the entire domain and time under statistically stationary conditions, thus yielding

(3.12) \begin{equation} 0 = -\nu \left \langle \overline {\frac {\partial u_i}{\partial x_{\!j}} \frac {\partial u_i}{\partial x_{\!j}}} \right \rangle _{V} + g \beta \left \langle \overline {u_2 T} \right \rangle _{V} + \left \langle \overline {u_i f_i} \right \rangle _{V} = - \left \langle \overline {\epsilon } \right \rangle _{V} + \left \langle \overline {P} \right \rangle _{V} + \left \langle \overline {u_i f_i} \right \rangle _{V} \end{equation}

where $\overline {\epsilon }, \overline {P} (=g\beta \overline {u_2 T})$ and $\overline { u_i f_i }$ indicate the time-averaged local energy dissipation rate, energy production by the buoyancy force and work done to the fluid by particles, respectively.

On the other hand, the kinetic energy equation of a particle can be obtained by multiplying $v^k_{p,i}$ using (2.5),

(3.13) \begin{align} \frac {\rm d}{{\rm d}t} \left ( \frac {v_{p,i}^{k} v_{p,i}^k}{2} \right ) & = v^k_{p,i} (F_{D,i}^{k} -g \delta _{i2}) \end{align}

(3.14) \begin{align} &= (v^k_{p,i} - u^k_i) {F}_{D,i}^{k} + u^k_{i} {F}_{D,i}^{k} - g v_{p,2}^{k} \end{align}

(3.15) \begin{align} &= -\, \frac {1}{\tau ^{\prime}_p} (u^k_i - v^k_{p,i})(u^k_i - v^k_{p,i}) + u^k_i {F}_{D,i}^k - g {v}_{p,2}^k , \end{align}

where the summation convention applies only to subscript $i$ . Equation (3.15) can be averaged over all particles in the domain and time under stationary conditions, yielding

(3.16)

with

(3.17) \begin{equation} \epsilon ^k_p = \frac {(u^k_i - v^k_{p,i})(u^k_i - v^k_{p,i})}{\tau ^{\prime}_p} , \end{equation}

where $\dot {n}_{in}$ is the number of particles inserted into the domain at each time step divided by $\Delta t$ . Here $(v^k_{p,2} )_{in}$ and $(v^k_{p,2} )_{out}$ are the vertical velocities of the particles inserted and eliminated at the top and bottom walls, respectively. Here $\epsilon ^k_p$ indicates the energy dissipation through friction between $k$ th particle and the fluid. Here $-g {v}_{p,2}^k$ is the energy produced by gravity acting on settling particles. The meaning of the tilde above each term (e.g. $\widetilde {v_{p,2}^{k}}$ ) is the average of all the particles in the domain and time. However, in $\widetilde {(v^k_{p,2} )_{in}}$ and $\widetilde {(v^k_{p,2} )_{out}}$ , the tilde only refers to the number of particles entering or exiting the channel. Here, $\tau ^{\prime}_p=\tau _p/(1+0.15 {\textit{Re}}_p^{0.687})$ and the variation in $\tau ^{\prime}_p$ over the particles is considered in the averaging process. The boundary conditions for the particles leaving the domain at the bottom and particles inserted into the domain at the top suggest that the left-hand side of (3.16) is zero. The balance equation for the particle energy (3.16) then becomes

(3.18) \begin{equation} 0 = - \widetilde {\epsilon ^k_p} + \widetilde {u^k_i F_{D,i}^k} + \widetilde {P_p}, \end{equation}

where $P_p =-g v_{p,2}^k$ is the production of particle energy owing to gravity.

From (2.4), the global energy transfer between the particles and fluid becomes

(3.19) \begin{equation} \left \langle \overline {u_i f_i} \right \rangle _V = - \varPhi _m \widetilde {u^k_i F_{D,i}^k}, \end{equation}

where $\varPhi _m (={N_p m_p}/{\rho V})$ corresponds to the mass loading of the particles, indicating that $\varPhi _m \widetilde {u^k_i F_{D,i}^k}$ is the energy loss owing to the work done by the particles. This indicates that $\left \langle u_i f_i \right \rangle _V$ is the rate of energy transfer from the particles to the fluids through a two-way interaction, and its sign is not predetermined. Table 4 lists the values of each term in (3.12) and (3.18) for $d_p/D=0.0015, 0.003$ and 0.0045, and $\varPhi _m=0.06, 0.09$ and 0.12, respectively. For the particle energy equation (3.18), each term is multiplied by $\varPhi _m$ for comparison with the terms in the fluid kinetic energy equation (3.12). The first noticeable difference between the kinetic energy equations for the fluid and particles is that the magnitude of the production or dissipation of particles is at least 300 and up to 3000 times larger than the terms of the fluid for the range of particle parameters in our problem. This is due to the large density ratio of $\rho _p/\rho =711$ . The balance, which is the sum of all three terms in each equation, ranges from 0.4 % to 8 % of the largest term for the kinetic energy equation of the fluid to 0.01 %–4 % for the kinetic energy equation of the particles. The relatively large balance for large particles was probably caused by a numerical error in the approximation of the delta function distribution of the feedback force in finite grids. The terms responsible for the energy transfer between the particles and the fluid, $\left \langle \overline {u_i f_i} \right \rangle _V$ in the fluid equation and $- \varPhi _m \widetilde {u^k_i F_{D,i}^k}$ in the particle equation, which are equivalent to each other, as suggested by (3.19), are comparable to each other, as highlighted in table except for the case with $d_p/D=0.0015$ . This mismatch is acceptable, given that the magnitude of $\varPhi _m \widetilde {u^k_i F_{D,i}^k}$ is extremely small and of the order of the balance. Another noticeable observation is that $\left \langle \overline {u_i f_i} \right \rangle _V$ is negative for $d_p/D=0.0015$ and positive for $d_p/D=0.003$ and 0.0045. This implies that the fluid is dragged by the settling large particles, whereas the local fluid motion at the small particle location is opposite to the average particle motion. Therefore, net energy is transferred to the fluid from large particles, whereas energy is transferred to small particles from the fluid on average.

Table 4. The value of each term in kinetic energy equations of fluid (3.12) and of particles (3.18) non-dimensionalized by $U_f^3/(\sqrt {Ra} D)$ at $Ra=10^6$ , where $\sqrt {Ra}(=1000)$ is introduced for better readability and easy comparison. The energy transfer between particles and fluid from both equations is highlighted for comparison. The percentage in the bracket denotes the magnitude of the balance relative to the largest term in the corresponding equation.

Figure 14. Horizontally averaged energy feedback defined in (3.19) non-dimensionalized by $U_f^3/D$ at $Ra=10^6$ . Coloured symbols are marked for the corresponding horizontal distribution presented in figure 15.

Figure 15. Energy feedbacks on the selected horizontal planes at $Ra=10^6$ : (a,d,g) at $y/D=0.05$ ; (b,e,h) at $y/D=0.5$ ; (c,f,i) at $y/D=0.95$ . $d_p/D=0.0015$ and $\varPhi _m=0.12$ for (a,b,c). $d_p/D=0.003$ and $\varPhi _m=0.12$ for (d,e,f). $d_p/D=0.0045$ and $\varPhi _m=0.06$ for (g,h,i). The energy feedback is averaged over a short time of $t U_f/D=0.1$ . Coloured symbols at the side of each panel are matched with the same symbols in figure 14 for convenience of comparison.

The detailed energy transfer between the particles and fluid can be investigated using the local information of $\overline {u_i f_i}$ . The distributions of the horizontally averaged energy transfer $\left \langle \overline {u_i f_i} \right \rangle _A$ for the two mass loadings $\varPhi _m=0.06$ and 0.12 is illustrated in figures 14(a) and 14(b), respectively. The energy transfer from the particle equation $- \varPhi _m \widetilde {u^k_i F_{D,i}^k}$ is also presented for comparison. Here, the wide tilde denotes the average of the particles at the same height. The oscillatory behaviour near the top wall was caused by the sudden introduction of a feedback force at the moment of release of a particle at the top wall, which was intermittent owing to the finite number of particles. As the particle size decreased, this phenomenon weakened, as the total number of particles increased to maintain the same mass loading. The average of these two energy transfer distributions over the gap height corresponds to both sides of (3.19) or the boldface numbers in table 4. The distributions of $\left \langle \overline {u_i f_i} \right \rangle _A$ and $- \varPhi _m \widetilde {u^k_i F_{D,i}^k}$ are similar for large particles with $d_p/D=0.003$ and 0.0045 and are positive everywhere. However, for small particles with $d_p/D=0.0015$ , the distributions of the two quantities show some differences near the top and bottom walls, where they are negative. The negative correlation between $u_i$ and $f_i (\simeq (u_i-v_{p,i})/\tau ^{\prime}_p)$ implies that $u_i$ and $v_{p,i}$ are opposite on average because $|v_{p,i}| \simeq v_t \gg |u_i|$ ; that is, the particle velocity is nearly equal to the settling velocity, which is much larger than the magnitude of the fluid velocity.

The distribution of the energy feedback on the selected horizontal planes, as shown in figure 15 can provide a more detailed energy transfer between the fluid and the particles. For the hexagonal cell structures for $d_p/D=0.0015$ and $\varPhi _m=0.12$ (figures 15 a, 15 b and 15 c), the negative $u_i f_i$ in the hot plume regions at the centre of each hexagonal cell and the positive $u_i f_i$ in the cold plume regions at the periphery of each cell at all heights are clearly due to the ascending and descending motions of the fluid in the corresponding regions, respectively. In particular, the local energy feedback at the midchannel (figure 15 b) is much stronger than the horizontally averaged value shown in figure 14(b), while the energy feedback near the bottom and top walls (figures 15 a and 15 c) is relatively weak, even resulting in a net negative value. Similarly, for the square cell structures with $d_p/D=0.003$ and $\varPhi _m=0.12$ (figures 15 d, 15 e and 15 f), negative $u_i f_i$ at the centre of each cell and positive $u_i f_i$ along the periphery are observed. The only difference is that not only in the core region but also in the near-wall regions, the net energy feedback is positive, as shown in figure 14(b). However, for the roll structures for $d_p/D=0.0045$ and $\varPhi _m=0.06$ (figures 15 g, 15 h and 15 i), positive $u_i f_i$ in the descending regions is more dominant than negative $u_i f_i$ in the ascending region at both the core and near-wall regions, although the vertical fluid motions in the ascending and descending regions are approximately symmetric. In summary, when the particles are large and heavy and hence the settling particles directly drag the fluid, more energy is transferred to the fluid from the particles in the descending region than the energy transferred to particles from the fluid in the ascending region. However, when the particles are small, the energy transfer between the settling particles and the fluid is comparable in both directions, resulting in a net energy transfer from the fluid to the particles.