1. Introduction
Thermal convection, which involves fluid motion induced by temperature gradients, is a ubiquitous and vital phenomenon in nature, with far-reaching implications for diverse fields of study, including fundamental sciences, technology and environmental flows. In the geophysical and astrophysical context, thermal convection plays a pivotal role in shaping the dynamics of the oceans, atmosphere and interior of celestial bodies such as stars and planets (Busse Reference Busse1978, Reference Busse1989; Wyngaard Reference Wyngaard1992; Mapes & Houze Reference Mapes and Houze1993; Maxworthy & Narimousa Reference Maxworthy and Narimousa1994; Atkinson & Wu Zhang Reference Atkinson and Wu Zhang1996; Getling Reference Getling1998; Marshall & Schott Reference Marshall and Schott1999; Thorpe Reference Thorpe2004; García-Melendo et al. Reference García-Melendo, Hueso, Sánchez-Lavega, Legarreta, del Río-Gaztelurrutia, Pérez-Hoyos and Sanz-Requena2013; Young & Read Reference Young and Read2017; Hanson et al. Reference Hanson, Duvall, Birch, Gizon and Sreenivasan2020; Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020). In oceanic flows, thermohaline convection drives the deep-ocean circulations (Marshall & Schott Reference Marshall and Schott1999; Rahmstorf Reference Rahmstorf2000), whereas in the atmosphere, an accurate estimate of thermally driven convection is crucial for weather predictions and climate calculations (Hartmann, Moy & Fu Reference Hartmann, Moy and Fu2001). Moreover, thermal convection operates in both the Earth's outer core and mantle (McKenzie, Roberts & Weiss Reference McKenzie, Roberts and Weiss1974; Cardin & Olson Reference Cardin and Olson1994; Christensen Reference Christensen1995; Zhong et al. Reference Zhong, Zuber, Moresi and Gurnis2000; Finlay & Amit Reference Finlay and Amit2011; Guervilly, Cardin & Schaeffer Reference Guervilly, Cardin and Schaeffer2019).
Emulsions (multiphase flows composed of two immiscible liquid phases with similar densities) play an important role in many contexts, from several industries to oil spills in oceans, e.g. the distribution of oil droplets is crucial for assessing environmental damage (Li & Garrett Reference Li and Garrett1998; French-McCay Reference French-McCay2004; Gopalan & Katz Reference Gopalan and Katz2010). Heat transfer is also important in these flows; therefore, herein we focus on exploring turbulent Rayleigh–Bénard (RB) convection in liquid–liquid emulsions using direct numerical simulations (DNS), which allow for a detailed analysis of fluid dynamics at the smallest scales. The RB convection, the buoyancy-driven flow arising from heating a fluid from below and cooling it from above, represents the most common configuration for turbulent thermal convection (Ahlers Reference Ahlers2009; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012; Shishkina Reference Shishkina2021). This has been extensively investigated for diverse geometries and various scenarios.
In thermally driven turbulent flows, such as turbulent RB convection, the thermal plumes serve as the primary carriers of heat (Ahlers Reference Ahlers2009; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012; Shishkina Reference Shishkina2021). Therefore, researchers have explored various approaches to enhance the overall transport of mass, momentum and heat within the flow by manipulating these coherent thermal plumes (Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012; Graham & Floryan Reference Graham and Floryan2021). Among several examples, we mention here surface roughness and grooved walls to facilitate the detachment of plumes from the boundary layers (Shen, Tong & Xia Reference Shen, Tong and Xia1996; Du & Tong Reference Du and Tong1998, Reference Du and Tong2000; Roche et al. Reference Roche, Castaing, Chabaud and Hébral2001; Qiu, Xia & Tong Reference Qiu, Xia and Tong2005; Stringano, Pascazio & Verzicco Reference Stringano, Pascazio and Verzicco2006; Shishkina & Wagner Reference Shishkina and Wagner2011; Tisserand et al. Reference Tisserand, Creyssels, Gasteuil, Pabiou, Gibert, Castaing and Chillà2011; Salort et al. Reference Salort, Liot, Rusaouen, Seychelles, Tisserand, Creyssels, Castaing and Chillà2014; Wei et al. Reference Wei, Chan, Ni, Zhao and Xia2014; Wagner & Shishkina Reference Wagner and Shishkina2015; Goluskin & Doering Reference Goluskin and Doering2016; Zhu et al. Reference Zhu, Stevens, Verzicco and Lohse2017; Jiang et al. Reference Jiang, Zhu, Mathai, Verzicco, Lohse and Sun2018; Zhang et al. Reference Zhang, Sun, Bao and Zhou2018), altering the wettability of the walls (Liu et al. Reference Liu, Chong, Ng, Verzicco and Lohse2022a), implementing a slippery surface (Huang et al. Reference Huang, Wang, Bao and He2022), applying geometry modifications (Huang et al. Reference Huang, Kaczorowski, Ni and Xia2013; Chong et al. Reference Chong, Yang, Huang, Zhong, Stevens, Verzicco, Lohse and Xia2017), the combination of inclination of the convection cell and confined geometries (Zwirner & Shishkina Reference Zwirner and Shishkina2018; Zwirner et al. Reference Zwirner, Khalilov, Kolesnichenko, Mamykin, Mandrykin, Pavlinov, Shestakov, Teimurazov, Frick and Shishkina2020) and the insertion of vertical partition walls into the convection cell (Bao et al. Reference Bao, Chen, Liu, She, Zhang and Zhou2015). Also, altering the fluid properties (Roche et al. Reference Roche, Castaing, Chabaud and Hébral2002; Silano, Sreenivasan & Verzicco Reference Silano, Sreenivasan and Verzicco2010), implementing pulsed heating power on the lower plate (Jin & Xia Reference Jin and Xia2008), adding polymer additives (Ahlers & Nikolaenko Reference Ahlers and Nikolaenko2010; Benzi, Ching & De Angelis Reference Benzi, Ching and De Angelis2010), incorporating shear (Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020; Wang, Zhou & Sun Reference Wang, Zhou and Sun2020) and rotating the convection cell (King et al. Reference King, Stellmach, Noir, Hansen and Aurnou2009; Zhong et al. Reference Zhong, Stevens, Clercx, Verzicco, Lohse and Ahlers2009b) have all been explored aiming to enhance heat transport. An alternative strategy to enhance the global heat transfer in thermal flows entails the introduction of a secondary phase. These encompass two immiscible fluids (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Lohse & Xia Reference Lohse and Xia2010; Chillà & Schumacher Reference Chillà and Schumacher2012), two immiscible fluid layers in the non-turbulent (Nataf, Moreno & Cardin Reference Nataf, Moreno and Cardin1988; Prakash & Koster Reference Prakash and Koster1994; Busse & Petry Reference Busse and Petry2009; Diwakar et al. Reference Diwakar, Tiwari, Das and Sundararajan2014) and turbulent regime (Xie & Xia Reference Xie and Xia2013; Yoshida & Hamano Reference Yoshida and Hamano2016; Liu et al. Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021, Reference Liu, Chong, Yang, Verzicco and Lohse2022b), liquid–liquid emulsions (Pelusi et al. Reference Pelusi, Sbragaglia, Benzi, Scagliarini, Bernaschi and Succi2021; Liu et al. Reference Liu, Chong, Ng, Verzicco and Lohse2022a), the generation of turbulence through the injection of air bubbles (Deckwer Reference Deckwer1980; Sekoguchi et al. Reference Sekoguchi, Nakazatomi, Sato and Tanaka1980; Sato, Sadatomi & Sekoguchi Reference Sato, Sadatomi and Sekoguchi1981; Tokuhiro & Lykoudis Reference Tokuhiro and Lykoudis1994; Deen & Kuipers Reference Deen and Kuipers2013; Dabiri & Tryggvason Reference Dabiri and Tryggvason2015; Gvozdić et al. Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018; Ng et al. Reference Ng, Spandan, Verzicco and Lohse2020) or through the formation of vapour bubbles via boiling (Oresta et al. Reference Oresta, Verzicco, Lohse and Prosperetti2009; Zhong, Funfschilling & Ahlers Reference Zhong, Funfschilling and Ahlers2009a; Schmidt et al. Reference Schmidt, Oresta, Toschi, Verzicco, Lohse and Prosperetti2011; Biferale et al. Reference Biferale, Perlekar, Sbragaglia and Toschi2012; Lakkaraju et al. Reference Lakkaraju, Stevens, Oresta, Verzicco, Lohse and Prosperetti2013; Guzman et al. Reference Guzman, Frkaczek, Reetz, Sun, Lohse and Ahlers2016a,Reference Guzman, Xie, Chen, Rivas, Sun, Lohse and Ahlersb; Wang, Mathai & Sun Reference Wang, Mathai and Sun2019), bubbles attached to the plate to mimic the boiling (Liu et al. Reference Liu, Chong, Yang, Verzicco and Lohse2022c), the inclusion of small particles (Oresta & Prosperetti Reference Oresta and Prosperetti2013; Park, O'Keefe & Richter Reference Park, O'Keefe and Richter2018), as well as non-colloidal suspensions experiencing both laminar (Kang, Yoshikawa & Mirbod Reference Kang, Yoshikawa and Mirbod2021) and turbulent flows (Chang & Ge Reference Chang and Ge2020; Demou et al. Reference Demou, Ardekani, Mirbod and Brandt2022).
In the case of emulsions, the breakup (or coalescence) of the fluid interface is expected to play a vital role in the heat transfer mechanism of RB convection flows. Previous studies have extensively investigated breakup and coalescence in turbulence (Deane & Stokes Reference Deane and Stokes2002; Villermaux Reference Villermaux2007; De Vita et al. Reference De Vita, Rosti, Caserta and Brandt2019; Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019b; Soligo, Roccon & Soldati Reference Soligo, Roccon and Soldati2019; Wang et al. Reference Wang, Mathai and Sun2019; Villermaux Reference Villermaux2020; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022). However, the presence of walls significantly impacts the behaviour of emulsion droplets, leading to the formation of clusters and complex structures (Scarbolo, Bianco & Soldati Reference Scarbolo, Bianco and Soldati2015). The classical Kolmogorov–Hinze theory (Kolmogorov Reference Kolmogorov1949; Hinze Reference Hinze1955) explains fluid breakup based on the balance between surface tension and inertial forces. However, in turbulent RB convection flows, buoyancy becomes a significant factor that can impact the breakup criteria (Liu et al. Reference Liu, Chong, Wang, Ng, Verzicco and Lohse2021), altering the heat transfer mechanism. This aspect remains an area of research that has not been comprehensively studied, and the present study aims to investigate it in detail. The objective of this work is, therefore, to fill this knowledge gap by comprehensively studying the behaviour of emulsions and their influence on heat transfer modulation in turbulent RB convection flows. We aim to gain a deeper understanding of emulsion dynamics near walls and their impact on heat transfer, which can provide valuable insights into the underlying mechanisms that govern the behaviour of emulsions in turbulent RB convection flows. These insights hold significance for optimizing and designing industrial processes that involve multiphase turbulent thermal convection.
In particular, we investigate the influence of key parameters, including the concentration of the dispersed phase and the dynamic viscosity and thermal diffusivity ratios between the two fluids, on the dynamics of the emulsions and their impact on overall heat transfer in the system. This research also addresses the knowledge gap on how interactions of dispersed droplets and their concentration affect the temperature and velocity fields of both carrier and dispersed phases, the internal energy budget and turbulent kinetic energy (TKE), leading to changes in global heat transfer. The study conducted by Demou et al. (Reference Demou, Ardekani, Mirbod and Brandt2022) focused on the impact of different volume fractions of rigid particles on heat transfer modulations in turbulent RB convection. Their findings showed that adding up to 25 % particle volume fraction to the single-phase flow enhanced the Nusselt number. However, surpassing this threshold led to a decrease in the heat transfer rate, attributed to the dense layering of particles near the wall. Building upon the insights of Demou et al. (Reference Demou, Ardekani, Mirbod and Brandt2022), we investigate a similar flow condition with a deformable secondary phase, represented by liquid–liquid emulsions instead of suspensions. This investigation provides a detailed analysis of the effects of the secondary phase deformability on heat transfer modulation, which, to the best of our knowledge, has not been previously explored. In this paper, since we did not employ any stabilizing mechanisms, the term ‘emulsion’ is not entirely accurate. Nevertheless, we use ‘emulsion’ together with ‘mixture’ for simplicity and consistency with previous literature (Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019; Yi, Toschi & Sun Reference Yi, Toschi and Sun2021; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022, Reference Crialesi-Esposito, Boffetta, Brandt, Chibbaro and Musacchio2023a).
The paper is structured as follows. Section 2 outlines the equations governing our problem, the employed numerical methods and a detailed description of our computational set-up with the list of our simulations. Section 3 features and analyses the significant findings of our work. This includes the Nusselt number variations, temperature and velocity statistics, a detailed description of heat transfer and TKE budgets, and the droplet-size distribution (DSD) analysis. In the last parts of § 3, we evaluate the effects of dynamic viscosity and thermal diffusivity ratios on DSD and heat transfer modulations. Conclusions and a future outlook are provided in § 4.
2. Methodology
2.1. Governing equations
To study emulsions in a turbulent RB convection, we introduce an indicator function 
$\mathcal {H}$, equal to 
$1$ in the volume 
$V_1$ occupied by the dispersed phase (fluid 
$1$) and 
$0$ in the volume 
$V_2$ occupied by the carrier phase (fluid 
$2$). The function 
$\mathcal {H}$ is governed by the transport equation
\begin{equation} \dfrac{\partial\mathcal{H}}{\partial\tilde{t}} + \tilde{\boldsymbol{u}}\boldsymbol{\cdot}\tilde{\boldsymbol{\nabla}}\mathcal{H} = 0, \end{equation}
where 
$\tilde {\boldsymbol {u}}=(\tilde {u},\tilde {v},\tilde {w})$ is the one-fluid velocity field, assumed continuous in the whole domain. Note that, henceforth, the symbol 
$\tilde {\boldsymbol {\cdot }}$ indicates a dimensionless scalar or vectorial quantity. The transport of 
$\mathcal {H}$ is coupled with the incompressibility constraint, the Navier–Stokes equations for a Newtonian fluid and the transport equation for the temperature field. These read, in dimensionless form, as
\begin{gather} \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}=0, \end{gather}
\begin{gather} \tilde{\rho}\left[\dfrac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{t}} + (\tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}}) \tilde{\boldsymbol{u}} \right] ={-}\tilde{\boldsymbol{\nabla}} \tilde{p}+\sqrt{\frac{Pr}{Ra}}\tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}\left[\tilde{\mu}\left(\boldsymbol{\nabla}\tilde{\boldsymbol{u}}+\boldsymbol{\nabla}\tilde{\boldsymbol{u}}^{\rm T}\right)\right] + \frac{\tilde{\boldsymbol{f}}_\sigma}{We} + \dfrac{\tilde{\hat{\rho}}\boldsymbol{e}_z}{Fr}, \end{gather}
\begin{gather} \tilde{\rho}\tilde{c}_p\left[\dfrac{\partial \tilde{\theta}}{\partial \tilde{t}} + (\tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}}) \tilde{\theta} \right] = \dfrac{\tilde{\boldsymbol{\nabla}}\boldsymbol{\cdot}(\tilde{k}\tilde{\boldsymbol{\nabla}}\tilde{\theta})}{\sqrt{{Pr}{Ra}}}. \end{gather}
In the above, 
$\tilde {p}$ is the hydrodynamic pressure, 
$\tilde {\theta }$ is the temperature and 
$\tilde {\boldsymbol {f}}_\sigma =\tilde {\kappa }\boldsymbol {n}_\varGamma \tilde {\delta }_\varGamma$ is the surface tension forces with 
$\tilde {\kappa }$ the interfacial curvature, 
$\boldsymbol {n}_\varGamma$ the normal vector and 
$\tilde {\delta }_\varGamma$ the Dirac-delta function (Scardovelli & Zaleski Reference Scardovelli and Zaleski1999). Note that 
$\boldsymbol {e}_z=(0,0,-1)$ is the unit normal vector oriented in the gravity direction. Equations (2.2), (2.3) and (2.4) are made dimensionless by introducing a reference length scale 
$L_s=H$ with 
$H$ the cavity height, the reference temperature difference 
$\Delta T=T_h-T_c$, i.e. the temperature difference between the bottom and top boundary, and a reference velocity taken as the free-fall velocity 
$U_f=\sqrt {\beta _r \mathrm {g} L_s \Delta T}$, where 
$\beta _r$ is the reference isothermal expansion coefficient and 
$\mathrm {g}$ is the module of the gravitational acceleration.
For each generic thermophysical property, we introduce a reference 
$\psi _r$, which is chosen in two different ways. In the first choice, 
$\psi _r$ is taken equal to the average property of the dispersed and continuous phase weighted by the total volume fraction 
$\varPhi$, i.e. 
$\psi _r = \psi _1\varPhi +\psi _2(1-\varPhi )$, where 
$\varPhi = (\int _V\mathcal {H}\,{\rm d}V)/V$ and 
$V=V_1+V_2$ the total volume of the domain. Following this definition, 
$\psi _r$ is also used to define the generic dimensionless thermophysical property 
$\tilde {\psi }$ (density 
$\tilde {\rho }$, dynamic viscosity 
$\tilde {\mu }$, thermal conductivity 
$\tilde {k}$ or specific heat capacity 
$\tilde {c}_p$) as 
$\tilde {\psi } = \psi /\psi _r$, where 
$\psi$ is computed with an arithmetic average, i.e. 
$\psi = \psi _1\mathcal {H}+\psi _2(1-\mathcal {H})$. Accordingly, 
$\tilde {\psi }$ can be finally expressed as
\begin{equation} \tilde{\psi} = \dfrac{\psi_1\mathcal{H}+\psi_2(1-\mathcal{H})}{\psi_1\varPhi+\psi_2(1-\varPhi)} = \dfrac{1+(\lambda_\psi-1)\mathcal{H}}{1+(\lambda_\psi-1)\varPhi}, \end{equation}
where 
$\lambda _\psi =\psi _1/\psi _2$ is the property ratio. In the second choice, 
$\psi _r$ is taken equal to 
$\psi _2$ and, therefore, 
$\tilde {\psi }$ in (2.5) simply reduces to 
$1+(\lambda _\psi -1)\mathcal {H}$.
Regardless of the employed approach to define 
$\psi _r$, the different dimensionless numbers in (2.3) and (2.4) are expressed as follows. First, we define the Prandtl number, 
$Pr=\nu _r/\alpha _r$, as the ratio of the reference viscous and thermal diffusivity and the Rayleigh number 
$Ra=\beta _r \mathrm {g} L_s^3 \Delta T/(\alpha _r\nu _r)$ to characterize the importance of buoyancy forces to the viscous forces. Next, we introduce the Weber number 
$We = \rho _rU_f^2L_s/{\sigma }$ as the ratio between the inertial and the surface tension forces with 
$\sigma$ representing the surface tension coefficient. Finally, we define the Froude number 
$Fr=U_f^2/(gL_s)$ as the ratio between the inertial and gravity forces. Note that 
$\hat {\tilde {\rho }}$ in the last term of (2.3) is the volumetric density field modified to account for the thermal effects in the gravity forces. By assuming that the flow is incompressible within the limits of the Oberbeck–Boussinesq approximation (Oberbeck Reference Oberbeck1879; Boussinesq Reference Boussinesq1903; Gray & Giorgini Reference Gray and Giorgini1976), 
$\hat {\rho }$ takes the following form:
\begin{equation} \tilde{\hat{\rho}} = \dfrac{1}{\rho_r}\left[\rho_{1} ( 1 - \beta_1 \Delta T \tilde{\theta} ) \mathcal{H} + \rho_{2} ( 1 - \beta_2 \Delta T \tilde{\theta}) (1-\mathcal{H})\right]. \end{equation}
An important dimensionless parameter is the Nusselt number 
$Nu$, i.e. the dimensionless heat flux, which serves as an indicator of the overall heat transfer rate within the RB cell. Here 
$Nu$ is defined as
\begin{equation} Nu = \frac{h L_s}{k}=\frac{\text{total heat flux}}{\text{conductive heat flux}}=\frac{-\displaystyle{\left[k \dfrac{{\rm d}T}{{\rm d}z}\right]}_{wall}L_s}{k_r \Delta T}, \end{equation}
where 
$h$ represents the convective heat transfer coefficient of the flow, 
$k$ is the local thermal conductivity of the emulsion and 
$k_r$ is the reference thermal conductivity taken equal to the average thermal conductivity of the emulsion.
Throughout this work, we adopt the first approach, where for each change in 
$\varPhi$ and 
$\lambda _\psi$ during the simulation campaign, 
$\psi _r$ is kept equal to the value for 
$\varPhi =0$, i.e. the single-phase configuration. This approach allows us to investigate turbulence and heat transfer modulation by fixing dimensionless parameters defined using the thermophysical properties of the entire emulsion, rather than those of one of the phases. For the sake of completeness and comparison with this first methodology, we perform two additional simulations using the second definition and report the results in the Appendix.
2.2. Numerical methodology
The governing equations (2.2)–(2.4) are solved on a uniform Cartesian grid with constant grid spacing in the three directions, 
$\Delta x = \Delta y = \Delta z$. The grid spacing is defined as 
$\Delta x = L_x/N_x$, 
$\Delta y = L_y/N_y$ and 
$\Delta z = L_z/N_z$, where 
$L_x$, 
$L_y$ and 
$L_z$ are the lengths of the computational domain and 
$N_x$, 
$N_y$ and 
$N_z$ the number of grid cells in the three directions. The so-called one-fluid formulation (Prosperetti & Tryggvason Reference Prosperetti and Tryggvason2009) is employed to discretize the governing equations so that only one set of equations valid for both phases is solved over the whole domain. The procedure is as follows.
First, the numerical algorithm defines a cell-averaged value of 
$\mathcal {H}$, which is called volume-of-fluid (VOF) function or volume fraction,
\begin{equation} \phi = \frac{1}{V_c} \int_{V_c}^{} {\mathcal{H}({\tilde{x}},\tilde{t})\,{\rm d}V_c}, \end{equation}
where 
$V_c=\Delta x\Delta y\Delta z$ is the volume of each computational cell. By applying the definition (2.8) to (2.1), the advection for 
$\phi$ reads
\begin{equation} \frac{\partial {\phi}}{\partial \tilde{t}} + \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} {(\tilde{\boldsymbol{u}}{\mathcal{H}})}={\phi} \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}. \end{equation}
The various VOF methods proposed in the literature differ in the way 
$\mathcal {H}$ is approximated. In the present work, the multi-dimensional tangent hyperbola interface capturing (MTHINC) algorithm developed by Ii et al. (Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012) is employed, in which 
$\mathcal {H}$ is approximated with a hyperbolic tangent as
\begin{equation} \mathcal{H}(\tilde{X},\tilde{Y},\tilde{Z}) \approx \tfrac{1}{2} +\left\{ 1+\tanh\left[\gamma_{th}(\mathcal{S}(\tilde{X},\tilde{Y},\tilde{Z}) + {d}_{th}) \right] \right\}, \end{equation}
where 
$(\tilde {X},\tilde {Y},\tilde {Z})$ is a local coordinate system, i.e. 
$\tilde {X}=(\tilde {x} - 0.5/L_s)/(\Delta \tilde {x})$, 
$\tilde {Y}=(\kern0.08em \tilde {y} - 0.5/L_s)/(\Delta \tilde {y})$ and 
$\tilde {Z}=(\tilde {z} - 0.5/L_s)/(\Delta \tilde {z})$ with 
$(\Delta \tilde {x}, \Delta \tilde {y}, \Delta \tilde {z})$ the grid spacing in the scaled directions of 
$\tilde {x}$, 
$\tilde {y}$ and 
$\tilde {z}$. In (2.10), 
$\gamma _{th}$ is a parameter controlling the sharpness of colour function (set equal to 
$2$ in the present work), 
$\mathcal {S}(\tilde {X},\tilde {Y},\tilde {Z})$ is the surface function and 
$d_{th}$ is the normalization parameter. The implementation details of the MTHINC method are reported in Ii et al. (Reference Ii, Sugiyama, Takeuchi, Takagi, Matsumoto and Xiao2012) and Crialesi-Esposito et al. (Reference Crialesi-Esposito, Scapin, Demou, Rosti, Costa, Spiga and Brandt2023c). Once the colour function is known, the local average thermophysical properties of the emulsion (e.g. density, dynamic viscosity, specific heat capacity and thermal conductivity) are updated using the local volume fraction, i.e. 
$\psi = \psi _1\phi + \psi _2(1-\phi )$.
Next, the momentum and the temperature equations (2.3), (2.4) are discretized on a regular Cartesian grid using a staggered arrangement, i.e. all the scalar fields are defined at the cell centres, except for the velocity components, which are defined at the corresponding cell faces (Harlow & Welch Reference Harlow and Welch1965). All the spatial derivatives are approximated with second-order central schemes and the equations are advanced with a second-order Adams–Bashforth scheme. The pressure equation is solved with a direct fast-Fourier-transform-based Poisson solver to impose exactly the incompressibility constraint (2.2). The simulations are performed using the open-source code FluTAS (fluid transport accelerated solver, https://github.com/Multiphysics-Flow-Solvers/FluTAS), which is parallelized using MPI/OpenMP directives in the CPU version and accelerated using OpenACC directives in the GPU version. FluTAS is capable of performing interface-resolved simulations of incompressible multiphase flows, optionally with heat transfer, as shown in several past studies (De Vita et al. Reference De Vita, Rosti, Caserta and Brandt2019; Rosti, De Vita & Brandt Reference Rosti, De Vita and Brandt2019a; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019b; De Vita et al. Reference De Vita, Rosti, Caserta and Brandt2020; Kozul et al. Reference Kozul, Costa, Dawson and Brandt2020; Scapin, Costa & Brandt Reference Scapin, Costa and Brandt2020; Cannon et al. Reference Cannon, Izbassarov, Tammisola, Brandt and Rosti2021; Dalla Barba et al. Reference Dalla Barba, Scapin, Demou, Rosti, Picano and Brandt2021; Rosti & Takagi Reference Rosti and Takagi2021; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022; Scapin et al. Reference Scapin, Dalla Barba, Lupo, Rosti, Duwig and Brandt2022; Crialesi-Esposito, Chibbaro & Brandt Reference Crialesi-Esposito, Chibbaro and Brandt2023b; Mirbod et al. Reference Mirbod, Abtahi, Bilondi, Rosti and Brandt2023; Scapin, Demou & Brandt Reference Scapin, Demou and Brandt2023). More details can be found in Crialesi-Esposito et al. (Reference Crialesi-Esposito, Scapin, Demou, Rosti, Costa, Spiga and Brandt2023c), where standard benchmarks and scaling tests for the CPU and GPU versions can be found.
2.3. Case description
Direct numerical simulations are carried out in the three-dimensional RB cell reported in figure 1. The cell is periodic along the two horizontal directions, with wall boundary conditions imposed at the bottom and top wall. The temperature of the top and bottom walls is uniform, constant and equal to 
$\tilde {\theta }_c = - 0.5$ and 
$\tilde {\theta }_h = 0.5$. Emulsions are enclosed between two infinitely long plates. The figure depicts a sketch of our computational domain with height 
$L_{\tilde {z}} = 1$ in the scaled wall-normal direction (
$\tilde {z}$) and plates’ dimensions 
$L_{\tilde {x}} = L_{\tilde {y}} = 2$ in the periodic streamwise and spanwise directions (
$\tilde {x}$ and 
$\tilde {y}$). In this work we want to study turbulent emulsions representative of an oil–water system at ambient temperature and in a normal gravity condition. Accordingly, we choose 
$Ra=10^8$, 
${Pr=4}$ and 
$Fr=1$. The Weber number 
$We$ is chosen equal to 
$6000$, which balances the requirements of avoiding excessive unphysical coalescence, which occurs at low 
$We$, and excessive breakup, which occurs at high 
$We$. Note that 
$We$ is a large-scale Weber number, and it is defined based on the height of the cavity and the free-fall velocity. As a preliminary analysis, we also computed a ‘small-scale’ Weber number based on the local droplet radius and local velocity, and we always found it of the order of 
$10^2$. Furthermore, the two phases share the same densities (
$\rho _1=\rho _2=\rho _r$), specific heat capacity (
$c_{p,1}=c_{p,2}=c_{p,r}$) and thermal expansion coefficient (
$\beta _1=\beta _2=\beta _r$). The viscosity ratio 
$\lambda _\mu =\mu _d/\mu _c$ and thermal diffusivity ratio 
$\lambda _\alpha =\alpha _d/\alpha _c$, however, vary in the range of 
$[0.1-10]$. To investigate the effects of the dispersed phase on heat transfer rate in RB convection, different volume fractions of the dispersed phase are considered, 
$0.0 \leq \varPhi \leq 0.5$. Table 1 summarizes the simulations performed in the present work.
Figure 1. Schematic of the three-dimensional turbulent RB convection cell with the Cartesian coordinate system. The domain dimensions along the 
$\tilde {x}$, 
$\tilde {y}$ and 
$\tilde {z}$ directions are (
$L_{\tilde {x}}$, 
$L_{\tilde {y}}$, 
$L_{\tilde {z}}$) 
$=(2$, 
$2$, 
$1)$. The liquid–liquid emulsions are heated from the bottom wall (depicted in red) and cooled from the top wall (in blue). The figure illustrates instantaneous isosurfaces of temperature and dispersed droplets a short time (
$\Delta \tilde t = 0.4$ time units) after adding droplets (
$\varPhi =0.2$) to the single-phase flow.
Table 1. List of simulations performed in this study.
To ensure that all spatial scales are appropriately resolved, the grid resolution is determined based on three criteria, as outlined in Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). First, the local mesh size must be smaller than the local Kolmogorov scale, 
$\eta _K(\boldsymbol {x},t)$, the local Batchelor scale, 
$\eta _B(\boldsymbol {x},t)$, and the local length scale of 
$\eta _T(\boldsymbol {x},t)$, where 
$\eta _K(\boldsymbol {x},t) = [\nu ^{3/4}\epsilon (\boldsymbol {x},t)^{-1/4}]$, 
$\eta _B(\boldsymbol {x},t) = \eta _K(\boldsymbol {x},t)Pr^{-1/2}$ and 
$\eta _T(\boldsymbol {x},t) = \eta _K(\boldsymbol {x},t)Pr^{-3/4}$. Here, 
$\epsilon (\boldsymbol {x},t)$ is the local kinetic energy dissipation rates per mass. Therefore, in order to meet this well-established criterion, the number of grid points in the wall-normal direction should satisfy 
${N_{\tilde {z}}} \ge {\epsilon _{max}(\boldsymbol {x},t)}^{1/4}(Pr/\nu )^{3/4}H$. Thus, for all of our cases, 
${N_{\tilde {z}}}^{min} = 440$. Additionally, the global mesh size should be smaller than the global length scales of Kolmogorov, Batchelor and 
$\eta _T$. To meet this requirement, the number of grid points in the wall-normal direction should satisfy 
${N_{\tilde {z}}} \ge [{Ra}{(Nu_{max}-1)}Pr]^{1/4}$. Thus, 
${N_{\tilde {z}}}^{min} = 352$, considering 
$Nu_{max}=39$ (the maximum 
$Nu$ achieved in our study). Overall, the minimum number of grid points in the wall-normal direction for our uniform grid is calculated as 
${N_{\tilde {z}}}^{min} = \max (440, 352) = 440$. Finally, in order to resolve all spatial scales inside the boundary layers, the resolutions within the thermal boundary layer (TBL) and velocity boundary layer (VBL) should meet the following conditions: 
${N_{\tilde {z}}}^{TBL} \ge \sqrt {2} (0.482) {Nu}^{1/2}{(0.982)}^{3/2}$ and 
${N_{\tilde {z}}}^{VBL} \ge \sqrt {2} (0.482) {Nu}^{1/2}{Pr}^{1/3} {(0.982)}^{1/2}$. This translates to a minimum requirement of 
${N_{\tilde {z}}}^{min, BL} = \max (5, 7) = 7$ grid points within the boundary layers. In this study, we employ a grid size of 
$1024\times 1024\times 512$ with uniform spacing in the 
$\tilde {x}$, 
$\tilde {y}$ and 
$\tilde {z}$ directions, ensuring a minimum of 10 grid points within the boundary layers. Consequently, all three resolution requirements mentioned above are satisfactorily met.
We first simulated the single-phase case (
$\varPhi = 0$), initiated from 
$\tilde {{\boldsymbol u}} = \tilde {\theta } = 0$, until a statistically stationary state. A random noise equal to 
$5\,\%$ of the prescribed temperature difference is superimposed on the initial temperature to promote a faster transition to the turbulent state. The multiphase simulations of different volume fractions 
$\varPhi$, denoted as cases 2–5, were initialized using the same initial velocity field, that is, the statistically steady solution of the single-phase case. Moreover, when investigating the impact of the viscosity ratio (cases 6–8), we started from the statistically stationary solution of cases 2 and 5 and modified the dynamic viscosity of both phases while maintaining the average dynamic viscosity of the emulsion constant. A similar procedure is followed for cases 9–11 where we examined the effect of the thermal diffusivity ratio. However, for the last two cases of 12 and 13, similar to cases 6–8, we examined the effects of viscosity ratio, but this time we did this by just varying the value of the dispersed-fluid viscosity while maintaining the carrier-fluid viscosity identical to that of the single-phase scenario. The results regarding cases 12 and 13 are presented in the Appendix.
Once a statistical steady state was reached, all simulations were continued under stationary conditions over a predefined time interval to collect the turbulent statistics. Specifically, at the statistical steady state, simulations were run at a fixed time step for an interval of approximately 
$\Delta \tilde {t}_{ss}\approx 500$. During the data collection period, we stored a substantial number of samples, 14 000 per case. Note that the time step used to advance the governing equations was always dynamically adjusted according to the Courant–Friedrichs–Lewy (CFL) condition with 
$\text {CFL} = 0.25$, except during the statistical sampling stage. Here, the time step was held constant and set equal to 90 % of the average time step at 
$\text {CFL} = 0.25$. To ensure the convergence of statistics, we systematically computed and averaged both first- and second-order statistics using varying sample sizes. We assessed the differences between these statistics over four distinct time intervals within 
${\Delta \tilde {t}_{ss}}$, i.e. (a) the first quarter, (b) the first half, (c) the first three quarters and (d) the entire duration of 
$\Delta \tilde {t}_{ss}$. The analysis revealed a progressive reduction in the differences between (c) and (d), rendering the difference negligible. For clarity, figure 2(a) illustrates this procedure for case 4, where droplets were added to the single-phase flow at 
$\tilde {t}=756$, and statistical sampling started at 
$\tilde {t}=1160$ once the simulation has reached the statistical stationary condition. In particular, the temporal evolution of the Nusselt number at the top and bottom walls is monitored starting at 
${\tilde {t}_{drop}}$, when droplets are added into a statistically stationary single-phase flow, and 
${\tilde {t}_{ss}}$, when the multiphase flow reaches a statistically stationary condition. The temporally averaged 
$Nu$ is computed over the time interval 
$[{\tilde {t}_{start}}, \tilde {t}]$, where 
$\tilde {t}\ge \tilde {t}_{start}$. Figure 2(b) presents the ratio between the time-averaged Nusselt number at the top and bottom walls, which consistently approaches unity once a statistically stationary condition is achieved.
Figure 2. (a) Temporal evolution (green and orange solid lines) and temporal average of the Nusselt number on the top and bottom walls for case 4; the time averaging starts from two different time instants: when adding the droplets, 
$\tilde {t}=756$ (blue and red solid lines) and 
$\tilde {t}=1160$ for the final statistical sampling (pink and dark-grey solid lines). (b) The ratio of temporally averaged Nusselt numbers of top and bottom walls.
3. Results and discussions
3.1. Flow statistics
In this section we present the findings of our investigation, starting with the case of emulsions with property ratios equal to 1 to focus on the modifications induced by the presence of an interface, only. Figure 3 illustrates the variation of the time-averaged Nusselt number, normalized by the corresponding value in a single-phase configuration, as a function of the droplet volume fraction. Our results indicate that the introduction of droplets into the single-phase flow enhances the Nusselt number, particularly at higher droplet volume fractions. The maximum enhancement observed is 
$10.2\,\%$ at 
$\varPhi =0.5$. To provide a more precise assessment of the Nusselt number variation, we fit a third-order polynomial to the simulation data, yielding the following polynomial function:
\begin{equation} \dfrac{\left\langle{Nu}_{t,b}\right\rangle}{\left\langle{Nu}_{sp}\right\rangle} = 1.0 + 0.07202\varPhi + 0.2598{\varPhi}^2 + 0.01632{\varPhi}^3. \end{equation}
Note that as 
$\varPhi$ approaches zero, 
$\langle {Nu}_{t,b}\rangle$ correctly converges to 
$\langle {Nu}_{sp}\rangle$. Furthermore, with increasing values of 
$\varPhi$, the linear, quadratic and cubic terms express an increase of the Nusselt number with 
$\varPhi$, indicating enhanced heat transfer due to the presence of the dispersed phase. It is interesting to note here that when deformable emulsions are replaced with rigid particles in the same RB convection flow, the Nusselt number exhibits a non-monotonic behaviour, as reported in the study by Demou et al. (Reference Demou, Ardekani, Mirbod and Brandt2022). In particular, the average heat transfer slightly increases by up to 
$\varPhi = 30\,\%$ and then decreases well below the single-phase reference value at 
$\varPhi =40\,\%$. This behaviour is attributed to the migration of particles toward the near-wall region, a mechanism absent in the case of emulsions. Here, we anticipate that the increase, 
$\langle {Nu}_{t,b}\rangle /\langle {Nu}_{sp}\rangle >1$, and the associated enhanced mixing, is due to the increase of small-scale turbulence induced by interfacial stresses, despite the decrease of the dispersed-phase concentration in the near-wall region. A comprehensive exploration of the turbulence modulation and of the changes of the diffusion and convection terms at various droplet volume fractions will be provided in detail in § 3.2.
Figure 3. Nusselt number at different volume fractions of the dispersed phase.
Figure 4 offers visual representations of instantaneous temperature isosurfaces at different dispersed-droplet volume fractions. The temperature fields qualitatively corroborate the results presented in figure 3. As the droplet volume fraction increases, we observe an increase in the thermal plumes originating from both the upper and lower plates, accompanied by a reduction of the size of the flow structures. To quantitatively confirm this conclusion, we conducted a quantitative assessment of thermal plumes, employing the established definition of thermal plumes condition provided by Liu et al. (Reference Liu, Chong, Ng, Verzicco and Lohse2022a) as
\begin{equation} |\theta - \langle \theta\rangle | > \sqrt{\langle(\theta - \langle\theta\rangle)^2\rangle}, \end{equation}
where 
$\theta (\tilde {x}, \tilde {t})$ being the local temperature and 
$\langle \ \rangle$ the spatial and temporal average. In figure 4(e) we presented the ratio of 
${\varSigma _{plume}}/{\varSigma _{plume, sp}}$, which indicates the volume fraction of the thermal plume for each case normalized by the volume fraction of thermal plumes in the single-phase case. Evidently, with an increased dispersed-fluid volume fraction, a noticeable rise (around 5 % at 
$\varPhi =0.5$) in thermal plumes from both the upper and lower plates is observed.
Figure 4. Instantaneous temperature isosurfaces at 
$\tilde {\theta } = \pm 0.1$ (blue colour corresponds to 
$\tilde {\theta } = - 0.1$ and red colour corresponds to 
$\tilde {\theta } = 0.1$) and for various dispersed-droplet volume fractions of (a) 
$\varPhi = 0$, (b) 
$\varPhi = 0.2$, (c) 
$\varPhi = 0.4$ and (d) 
$\varPhi = 0.5$. (e) Volume fraction of thermal plume for cases 1–5 (table 1) normalized by the volume fraction of thermal plumes in the single-phase case.
Figure 5(a–d) report the mean and root-mean-square (r.m.s.) temperature fields along the wall-normal direction for the carrier phase, the dispersed phase and the emulsion, where averages are taken in time and in wall-parallel planes. As one moves away from the hot (cold) wall, the average temperature gradually decreases (increases). This ultimately levels off to reach the temperature at the cavity centre. The temperature fluctuation profiles exhibit a clear pattern: starting from zero at the walls, the peak near the walls defines the edge of the TBL. Past this point, the fluctuations steadily decrease until they reach an almost constant value in the cavity's central region. To gain a deeper understanding of the heat transport mechanisms in turbulent RB convection, it is crucial to access the dynamics within the TBLs (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Chillà & Schumacher Reference Chillà and Schumacher2012; Li et al. Reference Li, Shi, du Puits, Resagk, Schumacher and Thess2012). In this region, the heat transport is mainly due to conduction. Hence, we calculate the dimensionless TBL thickness, denoted as 
$\tilde {\delta }_{t}$, as (Pope Reference Pope2000; Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Scheel & Schumacher Reference Scheel and Schumacher2014)
\begin{equation} \tilde{\delta}_{t}=\frac{{\delta}_{t}}{H} =\frac{1}{2} \left\langle{ \left |\frac{\partial \tilde{\theta}}{\partial \tilde{z}} \right| }_{wall}^{{-}1} \right\rangle= \frac{0.5}{Nu}, \end{equation}
which indicates an inverse relationship between TBL thickness and the temperature gradient at the walls. The TBL thickness is reported in figure 5(e) for all cases under investigation. In agreement with the data in figure 5(a), we note that adding droplets alters the wall-normal temperature gradient. Specifically, there is a noticeable increase in the temperature gradient near the wall as 
$\varPhi$ increases, which corresponds to a decrease in the TBL thickness. As also documented later (see figure 21b), the presence of small droplets in the near-wall region increases the local mixing and the global heat transfer. The data confirm thinner TBLs at a higher droplet volume fraction (
$\varPhi$), from 
$0.0143 \leq \tilde {z} \leq 0.0157$. At higher 
$\varPhi$ values, the flow experiences an increased level of heat transfer within the TBL, ultimately resulting in an improved total heat transfer rate. This conclusion can be confirmed by figure 6, which depicts the isosurfaces of thermal plumes (red colour) at 
$\tilde {\theta } = 0.3$ and dispersed droplets (grey colour) inside the boundary layer close to the hot wall, and for various dispersed-droplet volume fractions 
$\varPhi = [0\unicode{x2013}0.5]$. It is clear that adding the dispersed fluid to the single-phase flow enhances the mixing within the boundary layers and increases the amount of thermal plumes emitted from the plate and transported to the centre of the cavity, which finally improves the Nusselt number and total heat transfer. Additionally, it is observed that, at higher 
$\varPhi$, there are more regions where thermal plumes are enclosed by the dispersed droplets, which carry the thermal plumes to the central regions by convection.
Figure 5. (a,b) Mean 
$\langle \tilde {\theta }\rangle$ and r.m.s. 
${\tilde {\theta }}^{rms}$ temperature profiles of emulsion, (c,d) r.m.s. temperature profiles of carrier 
${\tilde {\theta }}_c^{rms}$ and dispersed 
${\tilde {\theta }}_d^{rms}$ phase, along the wall-normal direction for different droplet volume fractions; (e) TBL thickness as a function of droplet volume fraction. In (a) one of the subsets shows the mean temperature per phase for a region within the TBL.
Figure 6. A top view of the instantaneous temperature isosurfaces (red colour) in 
$\tilde {x}\unicode{x2013}\tilde {y}$ plane inside the hot-plate boundary layers at 
$\tilde {\theta } = 0.3$ and for various dispersed-droplet volume fractions of (a) 
$\varPhi = 0$, (b) 
$\varPhi = 0.2$, (c) 
$\varPhi = 0.4$ and (d) 
$\varPhi = 0.5$. Dispersed droplets are shown with a grey colour.
An alternative approach to determine the TBL thickness consists in finding the maximum r.m.s. of the temperature profiles, displayed in figure 5(b) for different droplet volume fractions. At higher 
$\varPhi$, the temperature fluctuations are more pronounced over the entire wall-normal direction, indicating increased mixing and heat transport at higher 
$\varPhi$ (as also observed in figure 6). Upon closer examination of the figure, it becomes evident that the point of maximum temperature r.m.s. slightly moves closer to the wall when increasing 
$\varPhi$, indicating a reduction in the TBL thickness, a trend consistent with the findings presented in figure 5(e). Moreover, figure 5(c,d) illustrate the temperature fluctuation profiles for each phase separately. Generally, higher fluctuations are observed at higher 
$\varPhi$ in the central regions of the cavity. An exception to this trend is observed within the carrier phase, where a slight reduction in fluctuation levels is discernible when 
$\varPhi$ exceeds 
$0.4$. Additionally, the dispersed-phase data (figure 5d) display a decrease of the temperature fluctuation peak at higher 
$\varPhi$, accompanied by a shift in the peak position away from the wall. To better understand these data, one crucial aspect deserving detailed analysis is the local distribution of the dispersed phase within the cavity. To compute this, we consider that the local volume fraction 
$\phi$ is defined as the portion of a computational cell occupied by the dispersed phase. This takes the value 1 when the cell is fully occupied by a droplet and zero when the cell is fully occupied by the carrier phase. The wall-normal distributions of the spatially and temporally averaged local volume fraction, denoted as 
$\langle \phi \rangle$, are displayed in figure 7 for different nominal values of the dispersed-phase volume fraction, 
$\varPhi$. An approximately uniform distribution is observed within the core of the cavity for all cases. Interestingly, however, fewer droplets can be found near the wall, except for the case of a binary mixture, where droplets are evenly distributed throughout the wall-normal direction as expected by symmetry considerations (
$\varPhi =0.5$ with density and viscosity ratio equal to 1 for the cases in the figure). Further information regarding the distributions and sizes of the dispersed droplets will be provided in § 3.4. In the case of turbulent RB of rigid-particle suspensions, Demou et al. (Reference Demou, Ardekani, Mirbod and Brandt2022) observed a distinct near-wall peak in the dispersed-phase distribution. This peak, associated with particle layering, becomes more pronounced at higher particle volume fractions. The authors attribute this layering to the strong wall-particle lubrication interaction, which stabilizes the wall-normal position of particles after reaching the wall. Consequently, it becomes increasingly challenging for particles within the first layer to disengage from it. The opposite trend is observed in the case of deformable droplets.
Figure 7. Mean wall-normal distribution of local droplet volume fraction for various droplet volume fractions.
Another global quantity of interest is the kinetic energy of the two phases. To investigate the contributions of various flow structures, it is useful to partition the average kinetic energy per unit mass into vertical and horizontal components as (Demou et al. Reference Demou, Ardekani, Mirbod and Brandt2022)
\begin{gather} \tilde{K}_c^h=\tfrac{1}{2} {\left[{({\tilde{u}}_c^{rms})}^2 + {({\tilde{v}}_c^{rms})}^2\right] }, \quad \tilde{K}_c^v=\tfrac{1}{2} {\left[{({\tilde{w}}_c^{rms})}^2 \right] }, \end{gather}
\begin{gather} \tilde{K}_d^h=\tfrac{1}{2} {\left[{({\tilde{u}}_d^{rms})}^2 + {({\tilde{v}}_d^{rms})}^2 \right] }, \quad \tilde{K}_d^v=\tfrac{1}{2} {\left[{({\tilde{w}}_d^{rms})}^2\right] }. \end{gather}
In particular, the vertical kinetic energy captures the dynamics of vertical motions within the core of the cavity, while the horizontal kinetic energy is linked to the velocity of the kinetic boundary layers (KBLs) near the walls. Figure 8 displays the variations in mean horizontal and vertical kinetic energy along the wall-normal direction for the various 
$\varPhi$ under investigation in both dispersed and carrier phases. As expected, the horizontal components (see panels a,b) exhibit a peak near the wall and gradually approach a nearly constant value in the cavity centre. This trend is similar to that of the temperature r.m.s. In contrast, the maximum vertical kinetic energy (depicted in panels c,d) is attained at the cavity centre, with a gradual decrease towards the wall. The location of the peak of the horizontal kinetic energy corresponds to the edge of the KBL. Note that, at higher 
$\varPhi$, the KBL remains relatively unchanged, whereas the TBL decreases. This suggests that the large-scale circulation structures do not significantly vary when 
$\varPhi$ changes. Given the value of the Prandtl number (
$Pr=4$ in this study), indicating a difference in velocity and TBLs, we expect the VBL to be thicker than the TBL. By comparing the locations of maximum horizontal kinetic energy and temperature r.m.s., our observations align with this expectation.
Considering the differences with the dispersed-phase volume fraction, figure 8(a,b) illustrate a noticeable damping of the horizontal components of the kinetic energy with the volume fraction 
$\varPhi$, which is observed in both phases, however, more evident for the carrier phase. Note, however, that horizontal velocity fluctuations do not directly contribute to turbulent heat transport. Conversely, as demonstrated in figure 8(c,d), the level of vertical fluctuations does not vary significantly with 
$\varPhi$, except for the binary flow at 
$\varPhi =0.5$, when we observe a reduction for both phases. Referring back to figure 7, we recall that droplets tend to be distributed predominantly within the central region of the cavity, with fewer near the cavity walls. Furthermore, by increasing the volume fraction of the dispersed droplets, there is an increase in their absolute concentration within the cavity core. This implies that at higher 
$\varPhi$, a larger number of dispersed droplets actively engages in the vertical large-scale circulations of the flow. At 
$\varPhi =0.5$, it is not possible to define the carrier and dispersed phase, and, indeed, the contribution to the vertical large-scale circulation becomes equal within statistical accuracy. When considering a suspension of rigid particles, Demou et al. (Reference Demou, Ardekani, Mirbod and Brandt2022) reported a significant decrease in both the horizontal and vertical components of both phases at higher volume fractions of the secondary phase, which indicates a weakening of the large-scale circulation structures. We have, therefore, seen that the increase in the Nusselt number is not associated with an increase in the large-scale motions (quantified by the average TKE). We therefore consider the energy spectra, in particular, the one-dimensional longitudinal spectra associated with the horizontal velocity (the same trends are observed for the other components). Results pertaining to the single-phase flow as well as emulsions with volume fraction 
$\varPhi = 0.2\unicode{x2013}0.5$ are displayed in figure 9 where we consider spectra 
$E_{hh}(\kappa _h)$ extracted inside the TBL and at the cavity midplan. As demonstrated, the presence of dispersed droplets notably enhances the energy at small scales (high wavenumbers), with this effect becoming more pronounced at higher 
$\varPhi$ values particularly noticeable within boundary layers (as depicted in figure 9a). This indicates higher 
$Nu$ at higher 
$\varPhi$ as 
$Nu$ is evaluated at the wall. Conversely, introducing the dispersed droplets reduces the energy at large scales (low wavenumbers), a phenomenon more distinctly observed at the centre of the cavity (figure 9b). As reported in previous studies, the presence of the interface provides an alternative mechanism for energy transfer at small scales, typically at scales smaller than the smallest of the corresponding single-phase flow (see Perlekar Reference Perlekar2019; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022). This enhances the small-scale mixing in the near-wall region and possibly explains the increase in global heat transfer in emulsions of two fluids with the same properties.
Figure 9. The time-averaged spectrum of the TKE as a function of the wavenumber at different droplet volume fractions (a) inside of TBL and (b) at the centre of the cavity. Wavenumbers are normalized by the lowest non-zero wavenumber 
$\kappa _0 = {\rm \pi}/ H$.
3.2. Heat transfer budgets
This section examines the heat transfer budget, which is derived by applying the phase-ensemble averaging to the heat transfer equation. Using the framework developed and employed in Marchioro, Tanksley & Prosperetti (Reference Marchioro, Tanksley and Prosperetti1999) and Zhang & Prosperetti (Reference Zhang and Prosperetti2010), the different contributions to the total heat transfer are written in terms of convective and diffusive fluxes in each phase (Ahlers et al. Reference Ahlers, Grossmann and Lohse2009; Ardekani et al. Reference Ardekani, Abouali, Picano and Brandt2018; Demou et al. Reference Demou, Ardekani, Mirbod and Brandt2022) as
\begin{equation} q_{tot}'' = C_c + C_d + D_c + D_d, \end{equation}where
\begin{gather} C_c ={-} (1-\langle \phi\rangle)\langle w_c'T_c'\rangle , \end{gather}
\begin{gather}C_d ={-}\langle \phi\rangle \langle w_d'T_d'\rangle, \end{gather}
\begin{gather}D_c = (1-\langle\phi\rangle)\alpha_{c}\left\langle\frac{{\rm d}T_c}{{\rm d}z}\right\rangle, \end{gather}
\begin{gather}D_d = \langle\phi\rangle\alpha_{d}\left\langle\frac{{\rm d}T_d}{{\rm d}z}\right\rangle. \end{gather}
Here, 
$C$ and 
$D$ refer to the convection and diffusion heat fluxes, respectively. Also, 
$w'=w-\langle w\rangle$ and 
$T'=T-\langle T\rangle$ are the wall-normal components of the dimensional velocity and temperature fluctuations. Given that 
$\alpha _d=\alpha _c=\alpha$ in cases 1–8 (refer to table 1), (3.6) and (3.7) can be re-expressed in the dimensionless form
\begin{equation} \tilde{q}_{tot}'' = \tilde{C_c} + \tilde{C_d} + \tilde{D_c} + \tilde{D_d}, \end{equation}where
\begin{gather} \tilde{C_c} ={-} \sqrt{Ra\,Pr}(1-\langle\phi\rangle)\langle\tilde{w}_c'\tilde{\theta}_c'\rangle, \end{gather}
\begin{gather}\tilde{C_d} ={-} \sqrt{RaPr}\langle\phi\rangle\langle\tilde{w}_d'\tilde{\theta}_d'\rangle, \end{gather}
\begin{gather}\tilde{D_c} = (1-\langle\phi\rangle)\left\langle\frac{{\rm d}\tilde{\theta}_c}{{\rm d}\tilde{z}}\right\rangle, \end{gather}
\begin{gather}\tilde{D_d} = \langle\phi\rangle\left\langle\frac{{\rm d}\tilde{\theta}_d}{{\rm d}\tilde{z}}\right\rangle. \end{gather}
The data of the heat transfer budget are reported in figure 10. In figure 10(a) we report the total convection and diffusion heat fluxes, confirming that diffusion dominates at the wall and approaches zero in the centre of the cavity, whereas convection is the mechanism active in the bulk. The sum of the two is constant in the wall-normal coordinate and equal to the total heat flux through the cavity, used in the definition of the Nusselt number, see figure 3, where the maximum enhancement was found to be about 10 % for the highest 
$\varPhi$ considered.
Figure 10. (a) The total convection and total diffusion heat transfer, and (b) the convection and diffusion heat fluxes for dispersed and carrier phases, along the wall-normal direction and for various droplet volume fractions. (c) The carrier- and dispersed-fluid convection heat fluxes at the centre of the cavity (
$\tilde {z}=0.5$), normalized by the total heat flux of each case with various 
$\varPhi$.
Figure 10(b) shows the variation of diffusion and convection heat fluxes for both carrier and dispersed phases along the wall-normal direction. As 
$\varPhi$ increases from 
$0$ to 
$0.5$, the relative contributions of carrier-phase convection and diffusion heat transfer decrease, whereas the dispersed-phase contributions steadily rise, eventually reaching the same as the carrier phase for the binary mixture, i.e. 
$\varPhi =0.5$.
To examine the contribution to the heat transfer of the carrier and dispersed phase independently of the increase in volume fraction, we display in figure 10(c) the convection heat transfer at the centre of the cavity (
$\tilde {z}=0.5$) for dispersed and carrier fluids, normalized by the total heat transfer of each case with different 
$\varPhi$. This shows that the contribution of the dispersed phase is not directly proportional to the volume fraction, 
$\varPhi$, but rather lower. This disparity is clarified by recalling that the dispersed-phase concentration is lower near the walls. Droplets tend to preferentially remain outside the TBL, resulting in lower temperatures than that of the hot wall (vice versa for the cold wall). Consequently, the relative contribution of the dispersed phase to the total heat transfer is smaller than that of the carrier phase, e.g. only 
$17\,\%$ for 
$\varPhi =30\,\%$.
3.3. Turbulent kinetic energy budgets
To quantify the energy production and transfer mechanisms in emulsions, we investigate the TKE budget. The derivation of the TKE budget starts from the transport equation for the fluctuating velocity 
$u_i'$ as
\begin{equation} \rho\left(\frac{\partial {u}_i'}{\partial t} + \frac{\partial {u}_i' {u}_j'}{\partial x_j}\right) ={-} \dfrac{\partial p}{\partial x_i} + \frac{\partial {\tau}_{ij}}{\partial x_j} + \sigma \kappa \delta_{\varGamma} + \rho \mathrm{g}\left[1 - \beta_{th} (T - T_0)\right], \end{equation}
where 
$\tau _{ij} = 2\mu D_{ij}$ denotes the stress tensor and 
$D_{ij}=({\partial u_i}/{\partial x_j} + {\partial u_j}/{\partial x_i})/2$ stands for the strain rate tensor. Upon multiplying (3.10) by the velocity fluctuation, 
${u}_i'$, and performing some algebraic manipulations, the TKE evolution equation can be expressed as
\begin{align}
\rho\left(\frac{\partial {u}_i' {u}_i'/2}{\partial t} +
\frac{\partial {u}_i' {u}_i' {u}_j'/2}{\partial x_j}\right)
&={-} \frac{\partial {u}_i' p}{\partial x_i} +
\frac{\partial {{u}_i' \tau}_{ij}}{\partial x_j} -
\tau_{ij}D_{ij}\nonumber\\ &\quad + \sigma \kappa
\delta_{\varGamma} {u}_i' + \rho \mathrm{g} [1 - \beta_{th}
(T - T_0)] {u}_i'.
\end{align}
Equation (3.11) can be handled in two ways: it can be time averaged over both phases, leading to the total TKE (3.12a), or it can be time- and phase-averaged with respect to a specific phase 
$m$ (either carrier or dispersed phase). This approach results in the TKE equation specific to one phase, reported in (3.12b) (Dodd & Ferrante Reference Dodd and Ferrante2016; Rosti et al. Reference Rosti, Ge, Jain, Dodd and Brandt2019b; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022). We have
\begin{gather} \frac{{\rm d}\mathcal{K}}{{\rm d}t} + \frac{{\rm d}\mathcal{T}_{j}}{{\rm d}\kern0.06em x_j} = \mathcal{P} - \varepsilon + {\varPsi}^{\sigma}, \end{gather}
\begin{gather}\frac{{\rm d}\mathcal{K}_m}{{\rm d}t} = \mathcal{P}_m - \varepsilon_m + \varXi_m^{v} + \varXi_m^{p}. \end{gather}
In (3.12a), 
$\mathcal {P}$ represents the rate of TKE production resulting from buoyancy forces, 
$\varepsilon$ denotes the dissipation of kinetic energy due to the viscous effects, 
${\varPsi }^{\sigma }$ is the power of the surface tension (due to droplet deformation, breakup and coalescence) and 
$\mathcal {T}_j$ is the flux responsible for the spatial redistribution of the TKE. Note that after volume averaging (3.11) over both phases, the transport terms vanish. Moreover, the second term on the left-hand side of (3.11) is omitted during phase averaging because, as discussed in Dodd & Ferrante (Reference Dodd and Ferrante2016) and Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019b), turbulent eddies cannot transport TKE across the interface of carrier and droplet fluids in immiscible fluids. Finally, the terms 
$\varXi _m^{v}$ and 
$\varXi _m^{p}$ in (3.12b) denote the viscous and pressure work rates on phase 
$m$, which represent the transport of TKE by viscous stresses and pressure, respectively. Note finally that under statistically steady state conditions 
${{\rm d}\mathcal {K}}/{{\rm d}t} = {{\rm d}\mathcal {K}_m}/{{\rm d}t} = 0$. To summarize, each term in (3.12a) and (3.12b) reads
\begin{gather} \mathcal{K} =\langle {u}_i' {u}_i'\rangle/2, \end{gather}
\begin{gather}\mathcal{P} = \langle\mathrm{g}[1-\beta(T-T_0)]{u}_i'\rangle, \end{gather}
\begin{gather}\varepsilon =\langle 2 \nu D_{ij}D_{ij}\rangle; \quad {\varPsi}^{\sigma} =\left\langle \frac{1}{\rho} \sigma \kappa \delta_{\varGamma} {u}_i'\right\rangle, \end{gather}
\begin{gather}\frac{{\rm d}\mathcal{T}_{i}}{{\rm d}\kern0.06em x_j} = \mathcal{P} - \varepsilon + {\varPsi}^{\sigma} = \frac{\partial \langle{u}_i' {u}_i' {u}_j'\rangle/2}{\partial x_j}, \end{gather}
\begin{gather}\mathcal{K}_m = {\langle{u}_i' {u}_i'\rangle}_m /2, \end{gather}
\begin{gather}\mathcal{P}_m = {\langle\mathrm{g}[1-\beta(T-T_0)]{u}_i'\rangle}_m, \end{gather}
\begin{gather}\varepsilon_m = {\langle 2 \nu D_{ij}D_{ij}\rangle}_m, \end{gather}
\begin{gather}\varXi_m^{v} = {\left\langle\frac{1}{\rho} \frac{\partial{u}_i' p}{\partial{x_i}}\right\rangle}_m, \end{gather}
\begin{gather}\varXi_m^{p} = {\left\langle\frac{1}{\rho} \frac{\partial{u}_i' \tau_{ij}}{\partial{x_j}}\right\rangle}_m. \end{gather}
Some considerations are needed on the surface tension term (
${\varPsi }^{\sigma }$). This accounts for the work done by the surface tension forces. It can be either a source of TKE (positive) or a sink of TKE (negative). It relates to the rate of change of the surface energy at the droplet interface and is inversely proportional to the rate of change of the droplet surface area, i.e. 
${\varPsi }^{\sigma }(t) \propto - ({1}/{We}) ({{\rm d}A(t)}/{{\rm d}t})$, meaning that an increase in droplet surface area leads to a decrease in TKE (
${\varPsi }^{\sigma }(t) < 0$) and vice versa (see also discussion in Dodd & Ferrante Reference Dodd and Ferrante2016). At steady state, when the surface area is on average constant, the work rate of the surface tension forces becomes zero for homogeneous isotropic and shear turbulence, as proved in Dodd & Ferrante (Reference Dodd and Ferrante2016) and Trefftz-Posada & Ferrante (Reference Trefftz-Posada and Ferrante2023). Here, one can replicate the derivation outlined in Trefftz-Posada & Ferrante (Reference Trefftz-Posada and Ferrante2023) and show that 
${\varPsi }^{\sigma } = 0$ holds for RB convection as well, given that the mean wall-normal velocity is zero at steady state. Finally, the integral across the cavity of the nonlinear transfer terms also goes to zero, as expected by the conservative nature of the nonlinearities of the Navier–Stokes equations.
We now discuss the results. Figure 11(a) presents the wall-normal distributions of TKE density at different volume fractions. As depicted in this figure, the TKE distribution exhibits a peak close to the edge of the boundary layers, associated to the horizontal velocity components, as shown in figure 8. Beyond this peak, we observe a gradual decrease with the approach of an approximately constant value within the central region of the cavity. Figure 11(b) reports bar charts of the volume-averaged TKE. As illustrated in both figure 11(a,b), introducing a dispersed phase dampens the TKE, resulting in a 20 % reduction at 
$\varPhi =0.5$, in agreement with the findings in Dodd & Ferrante (Reference Dodd and Ferrante2016) for isothermal emulsions in decaying turbulence. Despite the reduction, the overall heat transfer increases by about 10 %: this is attributed to the transfer of energy to small scales by the surface tension forces, which enhances near-wall mixing (cf. discussion of figure 9). Indeed, as shown in Crialesi-Esposito et al. (Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022, Reference Crialesi-Esposito, Chibbaro and Brandt2023b) for emulsions in homogeneous and isotropic turbulence, the presence of the dispersed phase provides an alternative mechanism for energy transfer at small scales and reduces the size of the smallest active flow scales. In terms of global heat transfer, this more than compensates for the reduction in TKE.
Figure 11. (a) Wall-normal profiles of average TKE and (b) volume-averaged bar charts of TKE at different dispersed-phase volume fractions.
According to figure 11(b), except for the case of the binary mixture, the kinetic energy density of the dispersed phase is lower than that of the carrier phase. This difference is due to the vertical component of the carrier-phase velocity (see figure 8) and related to the increased concentration of the carrier phase in the near-wall region. In other words, the carrier fluid has higher/lower temperatures than the dispersed phase on average and is, therefore, subject to stronger buoyancy.
The wall-normal profiles of the different terms of the TKE budgets are displayed in figure 12(a) in log scale to highlight the near-wall dynamics. Before discussing the behaviour of each term, we note that the differences between emulsions with different volume fractions are relatively small in the bulk of the flow and the trends of the TKE budget closely mirror the observations made for suspensions of rigid particles (Demou et al. Reference Demou, Ardekani, Mirbod and Brandt2022).
Figure 12. (a) The TKE budget terms along the wall-normal direction; (b,c) bar charts representing the volume average of the different terms in the TKE budget for various droplet volume fractions.
The data show that the TKE dissipation rate is highest in proximity to the wall, increasing with the volume fraction 
$\varPhi$, within the VBL. It reduces towards the centre of the cavity and eventually reaches an approximately constant value. In contrast, the TKE production rate, zero at the wall, increases towards the core of the cavity where it reaches a plateau, with values increasing at higher values of 
$\varPhi$. We also note that the work of the interfacial stresses increases with the amount of dispersed phase. This term provides energy to the near-wall region, at the expense of the kinetic energy in the centre of the cavity. This suggests that the near-wall dynamics is characterized by the relaxation and coalescence of smaller droplets in a laminar-like flow, while breakup dominates the dynamics of bigger droplets in the bulk (see also average droplet sizes across the cavity in figure 21). The reader is referred to Crialesi-Esposito et al. (Reference Crialesi-Esposito, Chibbaro and Brandt2023b) for an analysis of the relation between scale-by-scale energy transfer and droplet dynamics.
Figure 12(a) also provides insights into the wall-normal distribution of 
$-({{\rm d}\mathcal {T}}/{{\rm d}z})$, which represents the spatial redistribution term. As can be seen in the figure, nonlinear interactions transfer TKE from the central regions of the cavity (where 
$-({{\rm d}\mathcal {T}}/{{\rm d}z})<0$) towards the viscous regions near the wall (where 
$-({{\rm d}\mathcal {T}}/{{\rm d}z})>0$), where this is dissipated by viscous forces.
Next, we display in figure 12(b,c) the temporally averaged and volume-averaged TKE production and dissipation rates. The bar charts reveal that the TKE production due to buoyancy forces (panel b) is larger in the carrier phase, with values increasing with the volume fraction 
$\varPhi$ for both phases (except the special case of the binary mixture). This is explained by the increased buoyancy of the carrier phase, on average located closer to the walls, as discussed, explaining the larger vertical component of the kinetic energy. As mentioned above, the carrier phase is most likely to stay within the TBLs; its temperature approaches that of the nearby walls (see inset in figure 5a), and its density is, therefore, more likely to reach low/high values. In other words, the fastest rising plumes are expected to contain more of the carrier than of the dispersed phase, as suggested by the fact that the temperature–velocity fluctuations are larger in the carrier phase. To conclude, we also note an overall increase in the production with the volume fraction of the dispersed phase, of the order of about 10 %.
The viscous dissipation density is found to be only slightly larger in the carrier phase, indicating the presence of larger velocity gradients. Despite production being more pronounced in the carrier phase, the dissipation appears to be more equally distributed. Regarding the dissipation rate enhancement with 
$\varPhi$, as noted in Dodd & Ferrante (Reference Dodd and Ferrante2016), the addition of more droplets to the single-phase flow leads to an increase in the velocity gradient (
${\partial u_i}/{\partial x_j}$) near the droplet surfaces, resulting in higher viscous dissipation rates. This is in local balance with the work rate of the surface tension, as shown in Crialesi-Esposito et al. (Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022). Finally, for the case of a binary mixture, production and dissipation density are the same in the two phases, at least within statistical error, with the dissipation rate showing a slightly better statistical convergence.
The interphase energy transfer due to viscous and pressure forces are displayed in figure 13. As shown in Trefftz-Posada & Ferrante (Reference Trefftz-Posada and Ferrante2023), at steady state, when the work of the surface tension forces is on average zero, the weighted sum of the interfacial fluxes is zero,
\begin{equation} 0= (1-\phi) (\varXi_c^{v}+\varXi_c^{p}) + \phi( \varXi_d^{v}+\varXi_d^{p}). \end{equation}
The data in the figure show that viscous stresses transfer energy from the carrier to the dispersed phase, whereas the work by the pressure forces increases the TKE of the carrier phase at the expense of the dispersed one. For the case of a binary mixture, 
$\varPhi =0.5$, the different contributions vanish (within statistical error) by symmetry. Note, finally, that the viscous and pressure transport terms are significantly smaller than production and dissipation, as in the case of isothermal turbulence reported in Rosti et al. (Reference Rosti, Ge, Jain, Dodd and Brandt2019b) and Trefftz-Posada & Ferrante (Reference Trefftz-Posada and Ferrante2023).
Figure 13. Bar charts denoting the volume average of the (a) viscous and (b) pressure power terms in the TKE budget at various dispersed-droplet volume fractions.
3.4. Droplet-size distribution
The DSD measures the range of sizes that droplets can assume in a given multiphase system. Thermal convection may lead to a different distribution than homogeneous isotropic turbulence, as the interplay between flow patterns and dispersed droplets can alter both their sizes and overall distribution. In figure 14 we present visualizations of the instantaneous spatial distributions of dispersed droplets during the statistical stationary state for two different dispersed-droplet volume fractions, namely 
$\varPhi =0.2$ and 
$\varPhi =0.5$. Following the release of the dispersed phase into the system, droplets undergo recurring coalescence and breakup events. Eventually, when the system reaches a stationary state, it appears as depicted in figure 14, where at higher droplet volume fractions, dispersed droplets are more prone to coalescence, leading to the presence of larger droplets. This qualitative observation is further confirmed by examining the probability density function (p.d.f.) of droplet sizes.
Figure 14. Instantaneous distributions of dispersed droplets at (a) 
$\varPhi =0.2$ and (b) 
$\varPhi =0.5$. Dispersed droplets are coloured based on their temperature. For a clear visualization of the thermal plumes between the two plates and the corresponding droplets’ transport, the temperature range in the colourbar is restricted from 
$-0.1$ to 
$0.1$, i.e. any droplet with a temperature equal to or smaller than 
$-0.1$ and equal to or greater than 
$0.1$ is rendered with a uniform blue and red colour, respectively.
Figure 15(a) illustrates the p.d.f. of droplets of equivalent droplet diameter 
$d = (6V / {\rm \pi})^{1/3}$, where 
$V$ represents the droplet volume measured in the simulation, for different volume fractions, 
$\varPhi =[0.2-0.5]$. The equivalent droplet diameter is scaled by the cavity height. The analysis of the DSD aims to determine if it adheres to two different scaling laws: one proposed by Deane & Stokes (Reference Deane and Stokes2002) for small-size droplets (
$d^{-3/2}$) and the power law 
$d^{-10/3}$ introduced by Garrett, Li & Farmer (Reference Garrett, Li and Farmer2000) for large droplets. Based on the Kolmogorov–Hinze criteria, there exists a critical diameter, 
$d_{Hinze}$, such that droplets with 
$d< d_{Hinze}$ mainly experience breakup, while those with diameters greater than the Hinze length scale, 
$d>d_{Hinze}$, predominantly undergo coalescence. Therefore, the Kolmogorov–Hinze scale roughly determines the transition between the two regimes mentioned above, and it is typically estimated as (Hinze Reference Hinze1955)
\begin{equation} d_{Hinze} = \left(\frac{We_{cr}}{2}\right)^{3/5}\left(\dfrac{\sigma}{\rho_c}\right)^{3/5}\varepsilon^{{-}2/5}, \end{equation}
where 
$\varepsilon$ is the turbulent dissipation rate that can be computed in the RB configuration as suggested in Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010),
\begin{equation} \varepsilon = \dfrac{\nu_c^3}{L_s^4}(Nu-1)Ra{Pr}^{{-}2}. \end{equation}
Equation (3.15) requires a value for the critical Weber number to compute 
$d_{Hinze}$. For this study, we take 
$We_{cr} = [0.5-1.0]$, as suggested in Hinze (Reference Hinze1955), and based on this choice, the Hinze length scale lies in the range indicated by the two vertical lines in figure 15(a).
Figure 15. (a) The p.d.f. of the DSD for different droplet volume fractions. The solid-black lines represent the two scaling of 
$d^{-3/2}$ from Deane & Stokes (Reference Deane and Stokes2002) and 
$d^{-10/3}$ from Garrett et al. (Reference Garrett, Li and Farmer2000); (b) temporal evolution of the number of droplets, denoted as 
$N_{droplet}$, within the domain and (c) the time-averaged number of droplets 
$\langle N\rangle _{droplet}$ along with its associated fluctuation range, for distinct scenarios characterized by dispersed-phase volume fractions.
The DSD data confirm both power laws for small and large droplets at 
$\varPhi =0.2$ and 
$0.3$, with the smallest approximated Hinze length scale (
$\tilde {d}_{Hinze} = 0.0303$) providing a good estimate of the transition between the two scaling laws. However, the 
$-10/3$ law becomes less apparent in more concentrated cases, i.e. 
$\varPhi =0.4\unicode{x2013}0.5$. It should be mentioned that, for computing the droplet equivalent diameter 
$\tilde {d}$, all droplets are assumed to be spherical. This assumption ceases to be valid in the most concentrated cases, where significant deviations from sphericity occur due to the formation of large filaments that follow the thermal plumes. Consequently, a deviation from the 
$-10/3$ law is expected and observed in the p.d.f. graphs of these cases. Furthermore, in all scenarios, a secondary peak emerges at high values of 
$\tilde {d}$ (Mukherjee et al. Reference Mukherjee, Safdari, Shardt, Kenjereš and Van den Akker2019). This secondary peak indicates the presence of a few larger connected regions (large filaments) within the periodic simulation domain. Notably, this secondary peak becomes progressively more pronounced with increasing values of 
$\varPhi$.
Figure 15(b,c) illustrates the temporal evolution of the number of droplets, 
$N_{droplet}$, and the averaged number of droplets 
$\langle N\rangle _{droplet}$ (with their fluctuation range) for the different cases under investigation. In figure 15(b) the variation of 
$N_{droplet}$ is reported during the statistically stationary state. The fluctuations in the number of droplets, 
$N_{droplet}$, indicate a competition between droplet breakup and coalescence. The large dips in these graphs denote the coalescence of a significant number of droplets, which results in the formation of large droplets. These large droplets, however, are unstable and tend to breakup into smaller ones, which is reflected by the subsequent increase in 
$N_{droplet}$. Both figure 15(b,c) indicate that at higher dispersed-phase volume fractions, fewer droplets are present in the flow. This is because droplets are more likely to coalesce, forming larger droplets in cases with higher 
$\varPhi$. This observation aligns with the p.d.f. in figure 15(a), where a more pronounced second peak establishes at larger 
$\varPhi$ due to enhanced coalescence.
3.5. Effects of different viscosity ratios on the heat transfer
We now consider the influence of the viscosity ratio on the flow turbulence and heat transfer; see cases 6–8 in table 1. Specifically, we consider viscosity ratios 
$\lambda _{\mu }=(\mu _d/\mu _c)= 0.1$ and 10, while the other dimensionless parameters remain the same as for cases 1–5. Among the various dispersed-droplet volume fractions, we focus on 
$\varPhi =0.2$ and the case of a binary mixture 
$0.5$. As mentioned before in § 2.1, we define various dimensionless numbers, i.e. 
$Ra$, 
$We$, 
$Pr$ and 
$Fr$, based on an effective reference viscosity, which can be defined based on either the volumetric averaged viscosity or the viscosity of the continuous phase (traditional method). As a comparison between the two methods, it should be mentioned that, in defining dimensionless parameters, the traditional approach overlooks the thermophysical properties of the dispersed phase, masking certain effects, particularly at high values of 
$\varPhi$. For instance, increasing the viscosity of the dispersed phase tenfold simultaneously alters both the Rayleigh number in the dispersed phase (
$Ra_d$) and the viscosity ratio (
$\lambda _{\mu }$), obscuring the specific impact of viscosity ratio changes. The mixture rule, however, defines dimensionless parameters based on the volumetric averaged properties, isolating the effect of the viscosity ratio. Therefore, in this section we employ this mixture rule (cases 6–8, table 1) to compute the effective viscosity, while the traditional method is used for cases 12 and 13 and the related results are presented in the Appendix. The effective viscosity is calculated as the volume average of the viscosities of the two phases:
\begin{equation} \mu_{eff}=\varPhi \mu_d + (1-\varPhi) \mu_c. \end{equation}
In other words, the volumetric averaged viscosity, defined by the arithmetic average in (3.17), remains the same as the viscosity in the previous flow cases, whereas the viscosity of each phase varies. As an example, for viscosity ratio 
$\lambda _{\mu }=10$ at 
$\varPhi =0.2$, the viscosity of the dispersed phase is 
$2.8\mu _{1}$ and that of the carrier phase is 
$0.28\mu _{1}$, where 
$\mu _{1}$ is the viscosity of the cases with 
$\lambda _{\mu }=1$. Furthermore, note that the cases with 
$\varPhi =0.5$ and 
$\lambda _{\mu }=0.1, 10$ are identical; we therefore report results for only one case. With our choice, for the three ratios investigated, we expect the Nusselt number to stay the same if this is only a function of the effective emulsion viscosity: our results show that this is the case only for a viscosity ratio of 
$\lambda _{\mu }=0.1$.
Regarding the choice of using the linear definition that we adopted here to define the effective viscosity, it should be mentioned that, instead of a linear definition, a more precise definition (e.g. using a relation obtained from the emulsion rheological curves) can more effectively take into account the effect of the dispersed phase in computing the effective viscosity. However, we estimated the effective viscosity for different cases using the work of De Vita et al. (Reference De Vita, Rosti, Caserta and Brandt2019), and it is clarified that 
$(\mu _{eff}/\mu _{sp})_{estimated}$ varies in the range of 
$0.55$ to 
$1.35$. Despite the work of De Vita et al. (Reference De Vita, Rosti, Caserta and Brandt2019) for shear flows, the effective viscosity obtained using that work is not substantially different from the linear definition we employed, although they are not identical.
Figure 16(a,b) displays the effect of different dynamic viscosity ratios on the Nusselt number and the TBL thickness. For both volume fractions under examination, the Nusselt number varies weakly for 
$\lambda _{\mu } = 0.1$, while it increases for 
$\lambda _{\mu }=10$ and 
$\varPhi =0.2$. In particular, for 
$\varPhi =0.2$, the Nusselt number is increased by 
$25\,\%$ compared with the single-phase flow for viscosity ratio 
$\lambda _{\mu }= 10$, while it is enhanced only by 
$2.5\,\%$ and 
$1\,\%$ for viscosity ratios of 
$1$ and 
$0.1$. Consistently, the same trend is observed for the reduction in boundary-layer thickness; see panel (b). These observations suggest that a considerable heat transfer enhancement (around 24 %) can be achieved with a more viscous dispersed phase at the same Rayleigh number. Figure 16(c) reports the mean temperature for the different cases under investigation. We note that the temperature of the dispersed phase is slightly negative just outside the boundary layer closer to the hot wall, an effect more pronounced for the case with 
$\lambda _{\mu }=1$. This indicates that the transport of the droplet by plumes of the carrier phase occurs at a rate faster than the time required for heat to diffuse.
Figure 16. (a) Nusselt number, (b) TBL thickness and (c) mean temperature profiles along the wall-normal direction for the different various viscosity ratios and droplet volume fractions investigated; see legend.
The increase of the global heat transfer at viscosity ratio 
$\lambda _{\mu }=10$ and 
$\varPhi =0.2$ is explained by the increased turbulence in the carrier phase. As shown in figure 17, the horizontal and vertical components of the TKE are largest for the cases with a more viscous dispersed phase and 
$\varPhi =0.2$. In this case, the increase is most pronounced for the less viscous carrier phase; nevertheless, a more turbulent carrier phase induces more intense fluctuations also inside the more viscous dispersed phase. As discussed for the case 
$\lambda _{\mu }=1$, the majority of the heat transport is associated with the carrier phase as this penetrates deeper into the TBL. When this is relatively less viscous, turbulent fluctuations increase, and so does the heat transport associated with the carrier phase. Concerning emulsions with 
$\lambda _{\mu }=0.1$, we note only a weak decrease in the energy of the carrier phase and an increase of the turbulence in the dispersed phase, mainly in the bulk when comparing with the case 
$\lambda _{\mu }=1$. For the cases of binary mixtures (dotted lines in the figure), we notice an increase in the horizontal component of the less viscous phase inside the boundary layer.
Figure 18 displays the TKE spectrum for the flow with 
$\varPhi =0.2$ and at different locations: inside the TBL and in the core of the cavity. The results confirm an increase of the small-scale energy in emulsions when compared with the single-phase flow, an increase higher at viscosity ratios 
$\lambda _{\mu }=10$ and 1. In the case 
$\lambda _{\mu }=10$, we also observe an energy increase at large scales, which we attribute to the reduced viscosity of the carrier phase. This energy enhancement at large scales is consistent with the results of figure 17. Hence, we attribute the largest increase of the Nusselt number for 
$\lambda _{\mu }=10$ to the combined effect of increased small-scale mixing due to the presence of an interface, as for the case 
$\lambda _{\mu }=1$, and of increased turbulence in the carrier phase due to its reduced viscosity. The wall-normal profiles of the temperature fluctuations are reported in figure 19 for the two volume fractions under consideration, where values for each phase are presented. The data at 
$\varPhi =0.2$ show that the level of carrier-phase fluctuation increases close to the wall when this is less viscous, while those of the dispersed phase increase further from the wall when 
$\lambda _{\mu }=10$. In other words, the near-wall activity increases close to the wall for the case 
$\lambda _{\mu }=10$, as shown also by the horizontal velocity fluctuations in figure 17(a). The results pertaining to the binary mixtures are presented in figure 19(c,d) for the carrier and dispersed phase (the notation of the carrier and dispersed phase is maintained here to be able to identify the more or less viscous phases). The data indicate that the temperature fluctuation peak of the more viscous case moves towards the wall.
Figure 18. The TKE spectra 
$E_{xx}(\kappa )$ as a function of wavenumber at a droplet volume fraction of 
$\varPhi =0.2$ and different viscosity ratios (a) inside of TBL and (b) at the centre of the cavity. Wavenumbers are normalized by the lowest non-zero wavenumber 
$\kappa _0 = {\rm \pi}/ H$.
Figure 19. (a–d) Carrier- and dispersed-phase r.m.s. temperature profiles along the wall-normal direction at various viscosity ratios and droplet volume fractions.
Figure 20(a) reports the local wall-normal distributions of the dispersed phase at different viscosity ratios for the cases with 
$\varPhi =0.2$ and 
$0.5$. It can be seen that the decrease of the dispersed phase close to the wall observed in the case 
$\lambda _{\mu }=1$ remains when the dispersed phase is less viscous, 
$\lambda _{\mu }=0.1$, while the distribution is more uniform across the cavity when 
$\lambda _{\mu }=10$. This is attributed to the increase in turbulence and thinning of the TBL in the case of the more viscous dispersed phase. Also, at 
$\lambda _{\mu } = 10$ and 
$\varPhi =0.2$, we observe droplet layering close to the wall, as shown by a negative peak in the local distribution, roughly at the boundary-layer edge. In the case of a binary mixture, we observe that the less viscous phase tends to leave the wall, whereas the more viscous phase is preferentially found close to it.
Figure 20. (a) The local wall-normal distributions of the dispersed phase. (b) Number of droplets at the steady state condition. The blue and red dotted lines indicate 
$\varPhi = 0.2$ and 
$\varPhi = 0.5$, respectively (error bar indicating their fluctuation ranges). (c) The p.d.f. of the DSD for the cases with different viscosity ratios and 
$\varPhi$.
Next, we present the results regarding the DSDs. The steady-state number of droplets (as well as their fluctuation ranges) is depicted in figure 20(b), whereas the p.d.f. of the DSD is shown in figure 20(c). For the moderately concentrated case of 
$\varPhi =0.2$, as 
$\lambda _{\mu }$ increases, the turbulence intensity increases and the breakup becomes more frequent, leading to a higher number of smaller dispersed droplets and fewer larger ones (see the p.d.f. graphs of figure 20c). To gain insight on the DSD across the cavity, we display the average number of dispersed droplets and the average diameter of droplets along the wall-normal direction in figure 21 for the case with 
$\varPhi = 0.2$: the data reveal that more droplets are found closer to the wall for the cases with 
$\lambda _{\mu }=1$ and 10. In the bulk, we find more and smaller droplets in the case of a more viscous dispersed phase (
$\lambda _{\mu }=10$), a consequence of the more intense turbulence. Fewer and bigger droplets are observed in the bulk when 
$\lambda _{\mu }=0.1$. Interestingly, droplets of similar size are observed close to the wall for all viscosity ratios under investigation. The data shown in the figure confirm the presence of smaller droplets close to the wall, responsible for increased mixing, as discussed in the previous section for the flow cases at 
$\lambda _{\mu }=1$. The trends in figure 21 can be explained by the increased turbulent activity at 
$\lambda _{\mu }=10$ as well as depletion of the TBL for 
$\lambda _{\mu }=0.1$.
Figure 21. (a) The average number of dispersed droplets and (b) the average diameter of droplets, along the wall-normal direction for the case with 
$\varPhi = 0.2$ and viscosity ratios of 
$0.1$, 
$1.0$ and 
$10$.
For the binary mixtures (
$\varPhi =0.5$), the number of droplets is minimum at 
$\lambda _{\mu }=1$ and increases when the two fluids have different viscosity, 
$\lambda _{\mu }=0.1, 10$. The increase in droplet is here attributed to the preferential accumulation of the more viscous phase close to the wall and to an increase of the horizontal components of the TKE. In other words, decreasing 
$\lambda _{\mu }$ from 
$1$ to 
$0.1$ enhances the breakup rates, resulting in a higher number of smaller dispersed droplets with the less viscous ones migrating from the near-wall region towards the core of the cavity (see figure 20a).
To deepen our analysis, we report in figure 22 the heat transfer budget for the cases with different viscosity ratios, and 
$\varPhi =0.2$ and 
$0.5$. In figure 22(a) the bar charts indicate the wall-normal integral of the diffusion and convection heat fluxes, normalized by the total heat flux for each case. Figure 22(b,c) further details the wall-normal distributions of these heat flux terms, specifically focusing on the near-wall region, 
$\tilde {z} = [0\unicode{x2013}0.1]$. This range is selected as it encompasses the region where most variations occur.
Figure 22. (a) Wall-normal integral of heat fluxes transferred by diffusion and convection, normalized by the total heat flux of each case; (b,c) wall-normal distributions of convection and diffusion terms of both phases, normalized by the volume-averaged total heat flux of the single-phase case (
$\varPhi =0.0$) for volume fraction 
$\varPhi =0.2$ and 
$0.5$.
Starting with the emulsion with volume fraction 
$\varPhi =0.2$, the different plots confirm that the increase in heat fluxes at 
$\lambda _{\mu } = 10$ is due to the increase of the transport in the carrier phase, which is now less viscous than for 
$\lambda _{\mu } = 1$. The turbulent transport from the dispersed phase is similar among the different cases, being slightly lower in the case of viscosity ratio 1 (see in particular panel 22b). For the case of a less viscous dispersed phase, we note a decrease in transport from the carrier phase, which is only partially compensated by the increase in fluxes associated to transport by the droplets.
To conclude, we have shown that, at fixed Rayleigh number and effective viscosity, the Nusselt number increases when the viscosity ratio increases from 
$1$ to 
$10$ for volume fraction 
$\varPhi =0.2$, i.e. when the carrier phase becomes less viscous at 
$\lambda _{\mu } = 10$. In this case, in addition to the increase of small-scale mixing discussed for the case 
$\lambda _{\mu } = 1$, the turbulent fluctuations increase, the TBL is thinner and the more viscous droplets tend to remain closer to the wall. Our results indicate that the viscosity ratio needs to be considered explicitly when determining empirical correlations for the Nusselt number when larger than 1 (
$1 < \lambda _{\mu } \le 10$), while the idea of an effective viscosity can be used when the dispersed phase is less viscous (
$0.1 \le \lambda _{\mu } < 1$). Note finally that the addition of a more viscous phase to the same carrier fluid would lead to an increase of the effective viscosity. In this case, the increase in Nusselt number shown here will be modulated by the decrease in the effective Rayleigh number, so that both parameters, 
$\lambda _{\mu }$ and 
$Ra$, need to be considered. This has been shown in the Appendix where the traditional method has been employed, and a reduction in Nusselt number is observed at higher 
$\lambda _{\mu }$ due to the damped turbulence level caused by the reduction in the effective Rayleigh number.
3.6. Effects of various thermal diffusivity ratios on heat transfer
In this section we examine the effects of varying thermal diffusivity ratios (
$\lambda _{\alpha } = \alpha _d/\alpha _c$) on the heat transfer mechanism inside the cavity. We conduct simulations for three configurations (cases 9–11 in table 1) within the range of 
$0.1 \leq \lambda _{\alpha } \leq 10$ by changing the thermal conductivity of both phases (
$k_d$ and 
$k_c$). This is done so that the average thermal diffusivity of the emulsion remains constant and equals the thermal diffusivity of cases 1–5.
Figure 23 presents the normalized Nusselt number as a function of thermal diffusivity ratio for the three cases with the dispersed-phase volume fractions of 
$\varPhi =0.2$ and 
$0.5$. The behaviour differs between these two volume fractions. In particular, at 
$\varPhi =0.2$, the Nusselt number increases at smaller 
$\lambda _{\alpha }$, whereas, for the case of a binary mixture, it is maximum for 
$\lambda _{\alpha }=1$ and decreases when the thermal diffusivities of the two phases are not equal (
$\lambda _{\alpha }=0.1, 10$).
Figure 23. Nusselt number as a function of the thermal diffusivity ratio for the emulsion with a volume fraction of the dispersed phase 
$\varPhi =0.2$ and for the binary mixture.
To understand the effect of the thermal diffusivity ratio, we first display the wall-normal distributions of the dispersed phase and the total number of droplets; see figure 24. As shown in figure 24(a), we observe an increase in the local volume fraction 
$\phi$ close to the wall for 
$\lambda _{\alpha }=0.1$ and 
$\varPhi =0.2$, and a corresponding decrease in the core of the cavity. Conversely, when 
$\varPhi =0.2$ and 
$\lambda _{\alpha }=10$, the dispersed phase is, on average, very seldom in contact with the walls. To explain this, note that, when 
$\lambda _{\alpha }=10$, the thermal diffusivity of the carrier fluid decreases by nearly 65 %, while that of the dispersed phase increases by a factor 
$3.5$ (this is to keep the nominal average thermal diffusivity the same as for the emulsion with 
$\lambda _{\alpha }=1$). Thus, the dispersed phase absorbs/releases heat faster and leaves sooner the near-wall region due to buoyancy. Given the reduced local concentration near the wall (figure 24a), the local average diffusivity is closer to that of the carrier phase in the near-wall region, so effectively less than for the case with 
$\lambda _{\alpha }=1$. This is confirmed by the data in panel (c) of the same figure where we report the local average thermal diffusivity versus the wall-normal distance. Here, we indeed note a significant decrease close to the wall for the case 
$\lambda _{\alpha }=10$, when the dispersed phase very seldom reaches the wall, and an increase of the effective diffusivity for the case 
$\lambda _{\alpha }=0.1$. This is however not as pronounced as the decrease observed for the case 
$\lambda _{\alpha }=10$, which explains why the Nusselt number is only slightly increased when reducing the thermal diffusivity of the dispersed phase (
$\lambda _{\alpha }=0.1$). We speculate that this asymmetry is possibly due to the ratio between the time scale of thermal diffusion and transport: although the heat transfer is faster, the onset of motion is limited by the viscosity, which is constant among the different cases.
Figure 24. (a) The local wall-normal distributions of dispersed droplets. (b) The averaged number of dispersed droplets at the steady state condition (together with their transient fluctuation ranges), for the cases with different thermal diffusivity ratios and various 
$\varPhi$. The blue and red dotted lines indicate 
$\varPhi = 0.2$ and 
$\varPhi = 0.5$, respectively. (c) The local wall-normal distribution of thermal diffusivity.
As mentioned above, owing to the increased diffusivity, the dispersed phase more quickly adjusts to the local temperature and plumes quickly form. This reduces the dispersed-phase residence time in the near-wall region. Since this is also the region with high shear, we observe a reduction in the number of droplets for the case 
$\lambda _{\alpha }=10$ (see panel b of the same figure). The reduced presence of the dispersed phase in the near-wall region causes reduced small-scale mixing due to the fewer smaller droplets forming inside the boundary layers, as shown by the energy spectra in figure 25, where we observe less energy when 
$\lambda _{\alpha }=10$. Again, almost no difference is observed between emulsions with 
$\lambda _{\alpha }=1$ and 0.1.
Figure 25. The TKE spectra as a function of wavenumber at a droplet volume fraction of 
$\varPhi =0.2$ and different thermal diffusivity ratios (a) inside of TBL and (b) at the centre of the cavity. Wavenumbers are normalized by the lowest non-zero wavenumber 
$\kappa _0 = {\rm \pi}/ H$.
In the case of a binary mixture, see figure 24(a), again the more conductive phase is less likely to be found close to the wall and the average near-wall thermal diffusivity is effectively reduced, which explains the reduction in global heat transfer. As for the case with 
$\varPhi =0.2$, the reduction of the Nusselt number is caused by a reduction of the local average conductivity (see figure 24c); for a binary mixture, however, the number of droplets is almost independent of 
$\lambda _{\alpha }$ (see figure 24b).
The intensity of the wall-normal velocity fluctuations is displayed in figure 26 for both phases. The figure shows that the fluctuations reduce for 
$\lambda _{\alpha }=10$ when the dispersed phase more quickly adjusts to the surrounding temperature and the plumes more easily lose their drive. We note also that the cases with 
$\lambda _{\alpha }=1$ and 0.1 display similar values in Nusselt number and wall-normal velocity fluctuations, despite the dispersed phase is not more likely to be found close to the wall. The weak decrease in near-wall thermal diffusivity appears to be compensated by the weak increase in the number of droplets.
Figure 26. (a–d) Wall-normal profiles of the carrier- and dispersed-phase r.m.s. wall-normal velocity for the different thermal diffusivity ratios and droplet volume fractions investigated in this study.
These changes in the droplet distribution patterns and the corresponding values of the local diffusivity lead to the variations in the dispersed-phase convection flux shown in figure 27. Note that, since 
$\alpha _d \neq \alpha _c \neq \alpha$ when 
$\lambda _{\alpha } \ne 1$, the diffusion terms in the budget (3.7) take the following dimensionless form:
\begin{gather} \tilde{D_c} = (1-\langle \phi\rangle) \frac{\alpha_c}{\alpha} \left\langle\frac{{\rm d}\tilde{\theta}_c}{{\rm d}\tilde{z}}\right\rangle, \end{gather}
\begin{gather}\tilde{D_d} = \langle \phi\rangle\frac{\alpha_d}{\alpha}\left\langle\frac{{\rm d}\tilde{\theta}_d}{{\rm d}\tilde{z}}\right\rangle. \end{gather}
Figure 27. (a) Wall-normal integral of heat fluxes transferred by diffusion and convection, normalized by the total heat flux in 
$\varPhi = 0.0$; (b,c) wall-normal distributions of convection and diffusion terms of both phases, normalized by the volume-averaged total heat flux in 
$\varPhi =0.0$ and for cases with 
$\varPhi =0.2$ and 
$0.5$.
Figure 27(a) depicts the wall-normal integral of the convection and diffusion heat fluxes, normalized by the total flux in the single-phase flow. Figure 27(b,c), instead, illustrates the distribution of these fluxes along the wall-normal direction. As in figure 22, we focus on the region close to the wall (
$\tilde {z}=[0\unicode{x2013}0.1]$) as it exhibits the most significant variations. The budget terms in the figure confirm the reduction in near-wall diffusion when 
$\lambda _{\alpha }=10$, because of the reduced concentration of the dispersed phase in the near-wall region and the local reduction of the average thermal diffusivity. When normalizing the different contributions with the total heat flux, we observe a reduction of the contribution of the dispersed phase despite being more conductive because of the depletion of the near-wall region.
4. Conclusions and outlook
We have presented the results of interface-resolved simulations of liquid–liquid emulsions within a turbulent RB convection flow. Our focus lies in examining the effects of three primary properties on heat transfer rates and flow modulations in two-phase systems: (a) dispersed-droplet volume fraction within the range of 
$0 \leq \varPhi \leq 0.5$, (b) dynamic viscosity ratio in the range 
$0.1 \leq \mu _d / \mu _c \leq 10$, and (c) thermal diffusivity ratio in the range 
$0.1 \leq \alpha _d / \alpha _c \leq 10$. The remaining dimensionless numbers are kept fixed and equal to 
$Ra=10^8$, 
$Pr=1$, 
$We=6000$ and 
$Fr=1$. In the following, the key findings are summarized.
For fluids characterized by the same thermophysical properties, adding dispersed deformable droplets to the single-phase flow enhances the heat transfer rate in the cavity with respect to the single-phase counterpart. Differently from the case of solid particles, the enhancement is monotonic with the dispersed-phase volume fraction, 
$\varPhi$, with a maximum enhancement rate of 
$10.2\,\%$ for the case of a binary mixture.
To explain this finding, we recall that the presence of an interface is known to provide an alternative mechanism for energy transfer to small scales, extending the range of active flow structures towards even smaller scales (see Perlekar Reference Perlekar2019; Crialesi-Esposito et al. Reference Crialesi-Esposito, Rosti, Chibbaro and Brandt2022, among others). Hence, despite the fact that we also report a reduction of the average TKE, the increase of energy at the smallest scales is found to be responsible for increased mixing and, therefore, at the origin of the increase in the Nusselt number reported here. Note also that the large-scale mixing is almost unaltered, as determined by the cell size and global temperature difference.
Our results also show that the dispersed phase is less likely to be found closer to the walls, whereas the carrier phase is most likely to remain within the TBLs, even if the two fluids have the same thermophysical properties. Hence, the fastest-rising plumes contain more of the carrier than of the dispersed phase. As a consequence, analysis of the TKE and heat transfer budgets reveals that the turbulence production and heat fluxes are mainly associated with velocity–temperature fluctuations in the carrier phase.
As mentioned above, the carrier phase is most likely to stay within the TBLs, its temperature approaches that of the nearby walls, and its density is therefore more likely to reach low/high values. In other words, the fastest rising plumes are expected to contain more of the carrier than of the dispersed phase, as suggested by the fact that the temperature–velocity fluctuations are larger in the carrier phase.
We also consider fluids with different properties, particularly dynamic viscosity and thermal diffusivity, the latter changed by altering the thermal conductivity. The change in thermophysical properties is achieved through two different approaches. The first approach keeps the effective properties constant while changing the volume fraction of the emulsion and the property ratio. This choice allows for maintaining a constant Rayleigh number based on the average mixture viscosity. The second approach keeps the properties of the carrier fluid constant, changing only the dispersed-fluid properties. Throughout this work, the first approach is mainly used, but two additional cases are performed using the second approach for a change in dynamic viscosity only and discussed in the Appendix. In the first case, unlike the second, we observe an increase of the Nusselt number of about 25 % at 
$\varPhi =0.2$ when the dispersed phase is more viscous (
$\lambda _\mu =10$). This increase is attributed to two concurrent effects. The increase of small-scale mixing is due to the presence of an interface (as for the case 
$\lambda _\mu =1$) and to an increased level of turbulence in the less viscous carrier phase. In the case 
$\lambda _\mu =10$, we also observe an almost uniform distribution of the two phases across the cavity and an increased number of droplets despite these being more viscous. Small differences are observed between the emulsions with 
$\lambda _\mu =1$ and 
$0.1$, when the near-wall region is characterized by a reduction in the concentration of the dispersed phase. In the traditional (second) approach, on the other hand, we observe a significant damping in the turbulence level, leading to a 
$14\,\%$ reduction in the value of the Nusselt number when increasing the viscosity ratio from 
$0.1$ to 
$10$.
To investigate the effects of thermal diffusivity, only the first approach is employed throughout the paper. Based on the results, for higher thermal diffusivity of the dispersed phase, we report a lower residence time near the walls, with the dispersed phase almost never reaching the walls. Due to the faster heat diffusion, buoyancy is soon active, and the near-wall region is filled with the fluid with lower diffusivity. This leads to a significant reduction of the Nusselt number, about half that of the corresponding single-phase configuration, which is therefore attributed to a reduction of the local average diffusivity in combination with a reduced number of droplets in the depleted near-wall region, i.e. reduced small-scale mixing.
Finally, we have examined the DSD for the different cases. The results confirm the two scaling laws: 
$d^{-3/2}$ from Deane & Stokes (Reference Deane and Stokes2002) and 
$d^{-10/3}$ from Garrett et al. (Reference Garrett, Li and Farmer2000) for small and large droplets (Hinze criteria), characterized by the dominance of coalescence or breakup, for cases up to 
$\varPhi =0.3$. At higher 
$\varPhi$, a slight deviation from the 
$-10/3$ scale is observed due to the enhanced coalescence and the large deviation from the spherical shape of the dispersed phase. A better agreement with the scaling law could be observed by accounting for the droplet morphology in the p.d.f. calculation.
The configurations examined in this study may serve as a framework for future studies. In addition to direct comparisons with experiments, given the role of the near-wall distribution of the dispersed phase documented here, a potential follow-up investigation could consider thermally patterned walls with distinct structures and shapes on both the top and bottom plates. Varying the wetting properties of these walls could offer valuable insights into the interplay between the wetting properties of the walls and the modulation of heat transfer in the context of multiphase RB convection and provide a strategy to control the system behaviour, as shown in Liu et al. (Reference Liu, Chong, Ng, Verzicco and Lohse2022a). The role of the DSD on the heat transfer should also deserve further attention, as shown by the importance of surfactant for the effective viscosity of emulsions (Yi et al. Reference Yi, Wang, van Vuren, Lohse, Risso, Toschi and Sun2022). Moreover, in the present study we employed a linear relation for calculating the effective reference properties of the emulsion in order to define the non-dimensional numbers. Future work should focus on deriving a nonlinear relation for the effective viscosity from the emulsion's rheological curves (under turbulent RB convective flow) and using this nonlinear relation to better estimate the effective viscosity.
Funding
P.M. and A.M.B. were supported partially by the National Science Foundation (award no. 1854376). This research used resources from the Argonne Leadership Computing Facility (ALCF), which is a DOE Office of Science User Facility supported under contract DE-AC02-06CH11357.
Declaration of interests
The authors report no conflict of interest.
Appendix. Evaluation of the selected reference thermophysical properties
In this section we present and analyse the DNS results from two additional simulations. These simulations were conducted using the properties of the carrier phase as the reference for 
$\psi _r$, i.e. 
$\psi _r=\psi _2$. Accordingly, the dimensionless groups 
$Ra$, 
$We$ and 
$Pr$ are defined based on 
$\rho _c$, 
$\mu _c$, 
$k_c$ and 
$c_{pc}$, with changes in the viscosity or thermal diffusivity ratios achieved by tuning the corresponding properties in the dispersed phase only. Note that hereinafter, we refer to this approach as the traditional one. The two additional cases correspond to cases 
$12$ and 
$13$ in table 1 (
$\varPhi =0.2$, and 
$\lambda _{\mu } = 0.1, 10$) and are compared with cases 
$6$ and 
$7$ in table 1.
Figure 28 illustrates the Nusselt number, TBL thickness and mean temperature profiles for the cases 2, 12 and 13. By comparing the results of these three cases, it is observed in panels (a) and (b) that the Nusselt number (TBL thickness) decreases (becomes thicker) when we introduce a more viscous dispersed fluid (
$Nu$ decreases around 
$14\,\%$ when changing 
$\lambda _{\mu }$ from 
$0.1$ to 
$10$), which is an opposite behaviour compared with the method employed in the paper, which is based on keeping effective viscosity constant. This discrepancy arises because, in the constant-effective-viscosity approach, increasing the viscosity of the dispersed fluid is accompanied by a reduction in the viscosity of the carrier fluid, thereby ultimately enhancing the turbulence intensity of the emulsion, as clearly shown in figures 17 and 19). Conversely, in the traditional method, where the viscosity of the carrier fluid remains unchanged, the introduction of a more viscous dispersed fluid mitigates fluctuations. This damping effect is primarily observed both within the TBL, as depicted in figure 29 and also throughout the entire domain, as shown in figure 30. This last picture illustrates the profiles of temperature r.m.s. and the kinetic energy (both horizontal and vertical) of both fluids, respectively. Furthermore, in figure 28(c) a higher mean temperature is observed within the TBL at higher 
$\lambda _{\mu }$ values, indicating lower temperature gradients in close proximity to the wall, corresponding to an increase in the thickness of the TBL.
Figure 28. (a) Nusselt number, (b) TBL thickness and (c) mean temperature profiles along the wall-normal direction for the different viscosity ratios and at a moderate droplet volume fraction of 
$\varPhi =0.2$. Labels shown in legends with 
$^*$ are the cases obtained based on the traditional way.
Figure 29. (a) Carrier- and (b) dispersed-phase r.m.s. temperature profiles along the wall-normal direction at various viscosity ratios and at a droplet volume fraction of 
$\varPhi =0.2$.
In conclusion, it is important to note that in the traditional method, the dimensionless numbers of the carrier-fluid properties remain unchanged when the viscosity ratio varies. However, the same dimensionless numbers of the dispersed-fluid properties do change, which this method overlooks. For instance, the dispersed-fluid Rayleigh number (
$Ra_d$) decreases with higher 
$\lambda _{\mu }$ values. Therefore, if we calculate an effective Rayleigh number using the constant-effective-viscosity approach and define 
$Ra_{eff}$ as done in Liu et al. (Reference Liu, Chong, Ng, Verzicco and Lohse2022a), i.e. 
$Ra_{eff} = \varPhi Ra_d + (1-\varPhi )Ra_c$, we get a lower effective Rayleigh number, 
$Ra_{eff} < 10^8$. Conversely, the approach employed in our work allows us to keep the same 
$Ra_{eff}$ and change the viscosity ratio only. As a result, it offers the possibility to develop correlations for the Nusselt number 
$Nu$ as a function of the volume fraction 
$\varPhi$ based on the only parameter that is changed (
$\lambda _{\mu }$ in our case).
Open access
