Non-monotonic transport mechanisms in vertical natural convection with dispersed light droplets

Abstract We present results on the effect of dispersed droplets in vertical natural convection (VC) using direct numerical simulations based on a two-way fully coupled Euler–Lagrange approach with a liquid phase and a dispersed droplets phase. For increasing thermal driving, characterised by the Rayleigh number, Ra, of the two analysed droplet volume fractions, $\alpha = 5\times 10^{-3}$ and $\alpha = 2\times 10^{-2}$, we find non-monotonic responses to the overall heat fluxes, characterised by the Nusselt number, Nu. The Nu number is larger when the droplets are thermally coupled to the liquid. However, Nu values remain close to the 1/4-laminar VC scaling, suggesting that the heat transport is still modulated by thermal boundary layers. Local analyses reveal the non-monotonic trends of local heat fluxes and wall-shear stresses: whilst regions of high heat fluxes are correlated to increased wall-shear stresses, the spatio-temporal distribution and magnitude of the increase are non-monotonic, implying that the overall heat transport is obscured by competing mechanisms. Most crucially, we find that the transport mechanisms inherently depend on the dominance of droplet driving to thermal driving that can quantified by (i) the bubblance parameter $b$, which measures the ratio of energy produced by the dispersed phase and the energy of the background turbulence, and (ii) ${\textit {Ra}}_d/{\textit {Ra}}$, where ${\textit {Ra}}_d$ is the droplet Rayleigh number, which we introduce in this paper. When $b \lesssim O(10^{-1})$ and ${\textit {Ra}}_d/{\textit {Ra}} \lesssim O(100)$, the Nu scaling is expected to recover to the VC scaling without droplets, and comparison with $b$ and ${\textit {Ra}}_d/{\textit {Ra}}$ from our data supports this notion.


Introduction
Bubbles are ubiquitous. Within a liquid, they can play an important role in the transport of mass and heat. Such complex interactions of bubbles and liquids can be found in various applications and process technologies, for example in cooling systems of power plants, metallurgical industries, catalytic reactions and in the mixing of chemicals (Brennen 2005;Balachandar & Eaton 2010;Mathai, Lohse & Sun 2020). One commonly studied class of bubble-liquid interaction is the bubble column (Mudde 2005), where liquid turbulence is generated and sustained by a rising swarm of bubbles. This form of turbulence is typically referred to as pseudo-turbulence (Lance & Bataille 1991;van Wijngaarden 1998;Mercado et al. 2010) or bubble-induced agitation (Risso 2018).
Various parameters can be controlled to modulate heat transport in a bubbly flow. For instance, one can use microbubbles to increase heat transport in the boundary layer (Kitagawa & Murai 2013) or by inclining the domain (Piedra et al. 2015). The fluid properties can also be varied. For example, Deen & Kuipers (2013) studied the effects of bubble deformability and found localised increase of heat fluxes when bubble coalescence prevails in the near-wall region, whereas Dabiri & Tryggvason (2015) showed that nearly spherical bubbles tend to aggregate at the walls, which in turn agitate the thermal boundary layers and result in higher heat transport than for the case with deformable bubbles. From these studies, one key observation that can be made is that heat transport enhancement has been largely linked to boundary layer effects, e.g. thinning of the thermal boundary layers or ejection of thermal plumes. On the other hand, a recent experimental campaign using a homogeneous bubble column found that the heat transport, characterised by the Nusselt number Nu, not only increases by up to 20 times, but also becomes insensitive to the thermal driving of the background flow, characterised by the Rayleigh number Ra (Gvozdić et al. 2018;Gvozdić et al. 2019). The Ra-insensitivity persists across a range of bubble volume fractions α between 5 × 10 −3 and 5 × 10 −2 , implying that bubble-induced liquid agitation overwhelmingly dominates the heat transport mechanism across the thermal boundary layers. Indeed, the multifold enhancement in Nu is consistent with engineering estimates in the design of bubble column gas-liquid reactors (Deckwer 1980). Is there, however, any link between bubbly flows that directly influence the boundary layers versus bubble column experiments? And if any, are the boundary layers uniformly affected by the dispersed phase when α > 0? In this paper, we ask the question of how other parameters, specifically the density ratio of the dispersed phase to liquid phase, influence heat transport. Inspired by the water column experiments in Gvozdić et al. (2018) and to make contact with recent studies, we selected a set-up of natural convection in a rectangular cell containing a dispersed phase consisting of freely rising and deformable light droplets.
The model set-up of the flow is thermal natural convection, in particular, a flow sustained by applying a temperature difference between two opposing walls. Classical examples of thermal natural convection include Rayleigh-Bénard convection (Ahlers, Grossmann & Lohse 2009), where the hot wall is at the bottom and the cold wall at the top, and horizontal convection (Hughes & Griffiths 2008;Shishkina, Grossmann & Lohse 2016), where heating and cooling is applied at the same horizontal level. When the flow is confined between a hot vertical wall and a cold vertical wall, gravity acts orthogonal to the heat flux and this set-up is referred to as vertical natural convection (VC). For confined VC, the bulk flow is quiescent (see mean profiles in figure 1a and visualisation in figure 1b) and at low Ra, the laminar-like boundary layers are expected to dominate heat and momentum transport (Shishkina 2016). This flow is unlike the unconfined, doubly periodic VC (Ng et al. 2015(Ng et al. , 2017 where a mean shear is present and determines heat transport in the bulk flow region (Ng et al. 2018). Hereinafter, we refer to the rectangular VC cell set-up as VC, for simplicity. When light droplets are introduced into VC, we ask two specific questions: (i) Do the heat and momentum transport statistics exhibit monotonicity for droplets (i.e. when the density of droplets are close to the density of the liquid)? (ii) How important is the role of thermal coupling between the droplets and the liquid?
To answer these questions, we perform direct numerical simulations (DNS) of VC with droplets where we have control over the density ratio and thermal coupling of the droplet phase to the liquid phase. The droplets are fully coupled to the liquid phase DNS using the immersed boundary method (IBM) and the interaction potential approach, both of which are versatile numerical methodologies to simulate fully coupled fluid flows with deformable interfaces (e.g. Meschini et al. 2018;Spandan, Verzicco & Lohse 2018b;Viola, Meschini & Verzicco 2020). Furthermore, IBM offers some computational advantages over existing numerical methods for multiphase flows (e.g. volume of fluid, level-set and front tracking), for instance, the underlying discretised grid is fixed and no sharp interfaces need to be resolved (Spandan et al. 2017). Recent advancements in the numerical methodology have allowed the use of sparser discretisations of the deformable interface relative to the underlying grid (Spandan et al. 2018a) without compromising numerical accuracy, further easing the computational requirements for large-scale multiphase flows. The disadvantage of IBM, however, is that droplet coalescence or splitting is hard to model and correspondingly in this paper we refrain from attempting to do so.
Our paper is organised as follows: in § 2, we first describe the flow set-up and numerical details for the fluid and dispersed phase. In § 3, the numerical results are examined in detail. By analysing the near-wall heat fluxes ( § 4) and wall-shear stresses ( § 5), we relate the droplet driving dynamics to changes in the near-wall statistics. In § 6, we discuss and compare the influence of our selected parameters and the experimental parameters as reported in Gvozdić et al. (2018). Finally, in § 7, we summarise our results and provide an outlook.

Flow set-up
Our reference set-up is the single-phase VC flow (figure 1b), which is a buoyancy driven flow confined between two differentially heated vertical walls and two adiabatic horizontal walls. This reference flow will be referred to as the liquid carrier phase. The flow is governed by mass conservation, balances of momentum and energy conservation which within the Boussinesq approximation read, , 3) and repeated indices imply summation. In (2.2), f i is the back-reaction forces of the dispersed phase on the fluid arising from the IBM. At an Eulerian point k, it is defined as where N l is the number of Lagrangian markers associated with the Eulerian point, c l is a scaling factor that enforces conservation of momentum when transferring forces back and forth between the Lagrangian and Eulerian locations, Φ l k is the transfer function containing the shape function coefficients for each Lagrangian marker (here, based on the moving least squares (MLS) approximation) and F l i is the desired volume force component at the Lagrangians l (refer to de Tullio & Pascazio (2016) and Spandan et al. (2017) for details). For single-phase VC, f i = 0. The thermal analogue to f i in (2.2) is q θ in (2.3), which we selectively enable or disable in the present study. We define ρ ref as the reference density, θ ref as the reference temperature, β is the thermal expansion coefficient of the fluid, ν the kinematic viscosity and κ the thermal diffusivity, all assumed to be independent of temperature. The unit length is defined as the distance between the heated plates, L z , and the streamwise and spanwise domain lengths are L x = 2.4L z and L y = 0.25L z , respectively. Hereinafter, all length scales are non-dimensionalised by L z . No-slip and no-penetration boundary conditions are imposed on the velocity at all four walls, whereas periodic boundary conditions are imposed in the y-direction. The left and right walls are imposed with temperatures hotter and cooler than the reference temperature θ ref ≡ There are eight non-dimensional governing parameters for the VC flow with droplets and we have selected three parameters to vary, namely the strength of thermal driving Ra, the dispersed phase volume fraction α and the strength of droplet driving Ra d (Ra is defined below and Ra d is defined in § 2.5). The remaining parameters are fixed and consists of the Prandtl number (Pr, defined below), the domain aspect ratio (L x /L z ), the Weber number (We, defined in § 2.1), density ratio of the dispersed phase to fluid phase (ρ) and the ratio of droplet diameter to unit length (D/L z ).
The governing parameters for VC are the Rayleigh and Prandtl numbers which are, respectively, defined as where ≡ θ h − θ c . The aspect ratio (L x /L z ) can also be an additional control parameter for confined thermal convection problems (van der Poel, Stevens & Lohse 2011 ;Zwirner & Shishkina 2018), but at present, we restrict our analyses to a fixed value. Our simulations cover the values of Ra = 1.3 × 10 8 -1.3 × 10 9 and for Pr = 7, corresponding to water.
The typical flow response is described by the Nusselt and Reynolds numbers, which quantify the dimensionless heat flux and degree of turbulence, respectively. In (2.5a), f w ≡ −κ( θ xy /dz)| w is the wall heat flux and (·)| w denotes the wall value. Here, U s is the 'wind'-based velocity scale for VC (Ng et al. 2015) and accordingly, we set U s ≡ u xy,max , which is the maximum mean vertical velocity. Here, we use the notations · xy and · yz to denote x y-averaged and yz-averaged quantities, respectively (time-averaging is implied). The associated fluctuating components are denoted by (·) xy and (·) yz , e.g. u xy = u − u xy . With the addition of the thermal forcing term, q θ , in (2.3), a different definition for Nu becomes necessary because (d θ xy /dz)| z=0 / =(d θ xy /dz)| z=L z and the mean temperature equation now obeys w θ xy − κd θ xy /dz − q θ xy z = const. To overcome this difficulty, we employ the dissipation rate-based definition for the Nusselt number where ε θ is the volume-averaged thermal dissipation due to turbulent fluctuations and · denotes time-and volume-averaged quantities. When q θ = 0, (2.6) equals to (2.5a). The definition in (2.6) is also a direct analogue to the drag reduction calculations for multiphase Taylor-Couette flows (e.g. Sugiyama, Calzavarini & Lohse 2008; Spandan, Verzicco & Lohse 2018b), making it convenient when comparing heat transport at matched Ra (discussed in § 3.5). Throughout this paper, we will use (2.6) when reporting values of Nu, unless defined otherwise. The droplets are fully resolved using IBM for deformable interfaces and the interaction potential approach (de Tullio & Pascazio 2016;Spandan et al. 2017Spandan et al. , 2018a. The simulations are also coupled in a so-called four-way manner, i.e. the simulation is capable of handling droplet-fluid forcing, fluid-droplet forcing, droplet-droplet collisions and droplet-wall collisions (see § 2.1 for details on collision detection and modelling). Our numerical methodology differs from point-particle-type simulations with heat transport (e.g. Oresta et al. 2009): since the droplets with diameter D (at the point of injection) are significantly larger than the turbulent Kolmogorov length scale η, we therefore fully resolve the inhomogeneous hydrodynamic forces acting at the droplet interface. To illustrate this point, we wish to stress that D/η ≈ 7-19 in our simulations. Here, η ≡ (ν 3 /ε) 1/4 , where ε ≡ ν (∂u i /∂ x j ) 2 is the volume-averaged turbulent kinetic energy dissipation rate. The key points of our IBM are detailed in § 2.1. This is followed by numerical validations ( § 2.2), a description of the Lagrangian governing equations ( § 2.3), a description on the model for thermally coupled droplets ( § 2.4) and, finally, the droplet Rayleigh number ( § 2.5).

Numerical details
The liquid phase is solved using DNS by a staggered second-order accurate finite difference scheme and marched in time using a fractional-step approach (Verzicco & Orlandi 1996). We employ equal grid spacings in the x-and y-directions, whereas the z-direction is stretched using a Chebychev type clustering. The selected resolutions are constrained by considerations of three issues: (i) the resolution at the corner flow regions; (ii) the resolution at the bulk flow; and (iii) minimum number of grid points per droplet diameter.
Concerning point (i), we based our estimate from the minimum resolution guidelines proposed for laminar-like thermal convection simulations (Shishkina et al. 2010). As a check, a coarser simulation with 20 % fewer grid points results in Nu values that are within 0.5 %, indicating good convergence for our resolution. For point (ii), we determined that max[Δx + i ≡ Δx i /δ ν ] ≈ 2.4 (details in table 1), where Δx i are the grid spacings in each ith-direction and δ ν ≡ ν/u τ is the viscous length scale based on the local shear velocity scale u τ ≡ [ν(∂ u y (x)/∂z)| w ] 1/2 . Here, · y denotes averaging in the y-direction and in time. Point (iii) is closely related to point (ii); although the bulk resolutions are coarse, they are carefully selected such that D/(max[Δx i ]) 28, comparable to other immersed boundary studies in turbulent flow with finite-size particles (Wang, Vanella & Balaras 2019). Other numerical strategies are certainly possible, such as employing uniform grid spacings (Lu, Fernández & Tryggvason 2005) or by eliminating walls in the simulations (Uhlmann & Chouippe 2017), however, these strategies are either limited by the Reynolds numbers, or can be computationally costly. The resolutions employed here are therefore a careful compromise for our values of droplet rise Reynolds numbers, where U d is the time-averaged vertical rise velocity of the droplet. Additional validation tests for our IBM are provided in § 2.2. Finally, to justify this point, we compare our IBM resolutions with the minimum resolution conditions for a flow over a rigid sphere (Johnson & Patel 1999). Given that our maximum droplet rise Reynolds number, max[Re d ] ≈ 220, for an equivalent sphere Reynolds number, the dimensionless boundary layer thickness at its stagnation point is δ sp /D ≈ 1.13Re −1/2 d ≈ 0.08 (Schlichting & Gersten 2000). Our simulation resolution assures that at least two grid points reside within the droplet boundary layer. It may be tempting to treat this grid resolution as inadequate, however, we emphasise that this estimate is not only based on the extreme boundary layer criterion at the stagnation point, it is also based on the maximum Re d value and largest grid spacing in our set-up. Our IBM resolution improves at lower Re d (i.e. for thicker droplet boundary layers) and for finer near-wall grid spacings.
Our IBM employs the fast MLS algorithm (Spandan et al. 2017(Spandan et al. , 2018a. Two volume fractions are simulated: α = 5 × 10 −3 and 2 × 10 −2 (see table 1). We also fixed the ratio of initial droplet diameter to unit length, D/L z = 0.08. Therefore, at any point of our simulations, there are 12 droplets for α = 5 × 10 −3 and 48 droplets for α = 2 × 10 −2 . The ratio of D/L z employed in our simulations is somewhat larger than laboratory experiments with bubbles (e.g. Gvozdić et al. 2018) where the ratio of bubble diameter to channel width ∼O(10 −2 ). Sizes of dispersed phase have been shown to play a role in influencing momentum and heat transport (Shen, Ceccio & Perlin 2006;Kitagawa & Murai 2013;Verschoof et al. 2016); however, our choice of D/L z is necessary in order to avoid terminally expensive resolutions in accordance to our immersed boundary criteria (iii) above. To ascertain whether the periodic (spanwise) domain size is sufficiently large to avoid self-interactions of droplets, it is instructive to estimate the extent of a typical droplet's wake. As an approximation, the length of a sphere wake for a comparable value of Re d ≈ O(100) of our simulations is L w /D ≈ 1 (cf. Clift, Grace & Weber (2005), Chapter 5). Given that L y /D ≈ 3.1 > 1, the spanwise domain size is sufficiently large and we assume that the droplets do not interact with their own wake. and (n x , n y , n z ) = (1200, 120, 480) for Ra 0.7 × 10 9 . Here, T s /(L z /U ) and T s / t d are the total simulation sampling interval in terms of the free-fall velocity and droplet rise times, respectively.
At the start of each simulation, droplets are spatially initialised as spheres in a 4 × 3 array (in the xz-plane at y = 0.5L y ) for α = 5 × 10 −3 and a 12 × 4 array for α = 2 × 10 −2 . The droplet initial velocities are prescribed using a simple constant acceleration equation as a function of their vertical height, with the assumption that the droplet velocities are zero at x = 0. Droplets do not cross or touch the horizontal boundaries; once a rising droplet's interface is within D/2 from the top of the domain, the droplet is simultaneously removed and reinjected randomly at the bottom of the domain at the height of D. An alternative procedure would be to impose a stationary and homogeneous flux of droplets, which may be more realistic and closer to laboratory experiments. Indeed, in a real flow, there can be many bubbles and the number of bubbles is statistically constant. However, it is infeasible to obtain this number numerically by brute force. Moreover, by imposing a constant flux, the time scale of injection becomes an additional control parameter, which we want to avoid in order to keep our problem simple. In short, our reinjection procedure is a precise choice by simulation design which mimics the real system without introducing an additional parameter. However, the trade-off for maintaining a constant droplet volume fraction and constant number of Lagrangian markers is that the rate at which the droplets are injected fluctuates in our simulations.
During the reinjection of the droplets, the new (spherical) droplets are different from the removed (deformed) droplets. So, strictly speaking, this approach is not conserving and the horizontal walls can be interpreted as not adiabatic (although this is numerically imposed, see § 2). Nevertheless, additional analyses of the average magnitudes of heat flux in the vertical direction in the middle of the cell (not included in this paper) indicate that the values are considerably smaller by at most O(10 −1 ) of the horizontal laminar flux, /L z . Therefore, heat transport is more important and predominant in the horizontal z-direction as compared to the vertical x-direction, as will be discussed later in § 4.
The initial transient statistics at the start of the simulations are rapidly washed away after two to three droplet flow-through cycles, T s / t d where t d is the time-averaged droplet rise time. The reason for this is because of the multiple droplet removal and reinjection routines. Therefore, before recording statistics, we conservatively discarded a minimum of five droplet flow-through cycles, which correspond to discarding the first 50 to 150 sampling intervals depending on Ra (defined by T s /(L z /U ), where U ≡ (gβΔL z ) 1/2 is the free-fall velocity). Thereafter, at least 20 droplet rise intervals are recorded for each simulation. The resulting averaged droplet spatial distributions in our simulations are uniform, as shown in figure 3(b). In total, approximately 2 M central processing unit (CPU) hours were consumed.
As introduced earlier, we refrain from simulating droplet coalescence and splitting since these phenomena are extremely challenging tasks from both a physical and numerical point of view. Instead, we use 2562 Lagrangian markers (equivalent to 5120 triangulated faces) to represent each discretised droplet interface. This approach is computationally efficient and scalable since it eliminates the need to stitch or regenerate meshes but is, however, an inherent limitation of the present IBM-interaction potential approach.
Collision events (such as when two droplets get close to each other) are exceedingly rare even from estimates of our α = 2 × 10 −2 cases. The reasons are because the droplets are randomly injected into the flow without any overlap, rise almost vertically, and are not strongly swept by the background large-scale circulation unlike other flow set-ups with a strong mean shear (e.g. Spandan et al. 2018b). Nevertheless, we still employed the collision detection algorithm of Spandan et al. (2018a) in our numerics and the elastic potential collision model by Spandan et al. (2018b). Briefly, when two or more Lagrangian markers from different droplets reside in the same Eulerian cell at any time step, the elastic potential repulsive force is applied to all Lagrangian markers in the Eulerian cell. This force is proportional to the square of the inverse distance between the marker and the centre of the Eulerian cell.
Upon reinjection, the initial interfacial temperature of the droplet, θ k init , is equated to the averaged temperature of the immersed fluid, according to where θ k is the fluid temperature value interpolated on the barycentre of each triangulated face using MLS and N f is the total number of triangulated faces forming the discretised droplet interface; θ k init is subsequently forced using IBM to the Eulerian grid. A small initial droplet vertical velocity ∼O(10 −2 )U is also prescribed.
For the droplet boundary conditions, we assume that the droplets have negligible thermal inertia and are surfactant-laden. The first assumption implies a small droplet Biot number, defined by Bi ≡ l d h/k d (where l d is the characteristic droplet length scale, h is the convective heat transfer coefficient and k d is the thermal conductivity of the droplet interface) so that the internal droplet temperature can be approximated by a uniform temperature in accordance with the lumped-capacitance model (Wang, Sierakowski & Prosperetti 2017). Our reasoning for the small droplet Biot number value is as follows: by substituting the length scale l d ≡ D/6 (ratio of sphere volume to sphere surface area) and the heat transfer Here, Nu d is the droplet Nusselt number and k f is the thermal conductivity of the fluid. Next, assuming small Péclet values and neglecting phase change, we approximate Nu d ∼ O(1) (cf. Oresta et al. 2009). Finally, we assume k d > k f , and so Bi < 1. The temperature boundary condition at the droplet interface is then a homogeneous time-dependent Dirichlet boundary condition (the thermal model is discussed in § 2.4). The second assumption implies that the droplet boundary conditions are no-slip and impermeable for velocity. Differences could be expected for free-slip boundaries corresponding to clean interfaces: free-slip boundaries prevent viscous boundary layers from forming and therefore would have negligible contributions to local viscous dissipation. However, since clean bubbles present a unique set of challenges to achieve in laboratory environments, we assume, for the sake of simplicity, the extreme scenario where the bubble interfaces are saturated with surfactants. Indeed, for physical systems with surface-active impurities, droplet interfacial dynamics may be closely approximated by a no-slip interface (Duineveld 1995;Jenny, Dušek & Bouchet 2004). These simplified boundary conditions also have the added benefit that they can be handled easily from a numerical point-of-view, and hence, are computationally efficient given the size of the flow problem.
Owing to deformation, individual droplet volumes can vary slightly throughout the simulation, but fluctuate about a constant reference volume -this is the underlying approach of the interaction potential (IP) model (described in § 2.3). To quantify the droplet deformability, we define the Weber number, We ≡ ρ ref U 2 D/σ , which yields the ratio of inertia to capillary forces, where σ is the surface tension. In our simulation strategy, σ is not prescribed explicitly. Rather, an additional tuning step is performed to obtain a set of IP model parameters such that We ≈ 3 × 10 −2 . The tuning step consists of the following: a set of model parameters is first estimated from existing simulations, for instance, from Spandan et al. (2018b). Then, with the selected model parameters, we performed controlled test simulations of one droplet interacting with simple flows for which reference analytical and computational data are available, specifically, a droplet in shear flow (e.g. Maffettone & Minale 1998) and a droplet in cross-flow (e.g. Loth 2008; Schwarz, Kempe & Fröhlich 2016). Finally, in accordance with the tuning criteria described in Spandan et al. (2017), σ is reverse-engineered by matching the droplet deformation dynamics to the reference results in Maffettone & Minale (1998) and Schwarz et al. (2016). It is emphasised that in order to simplify existing continuum models, this tuning step is a necessary and felicitous step in the implementation of our numerical model. We have also chosen to simulate a constant We value. The reason for this is because, for our selected parameters, the background buoyancy driven flow is stronger than droplet-induced agitations (discussed later in § 6). Therefore, we expect deformability to play a weaker role than other governing parameters for the droplets, for instanceρ, D/L z or α. After extensive precursor simulations and checks, we decided to simulate droplets at half the density of the fluid, i.e.ρ ≡ ρ d /ρ ref = 0.5, which is within the numerical stability limit for explicit IBM time integration schemes (Schwarz, Kempe & Fröhlich 2015). The value ofρ is kept constant throughout our simulations. Another reason why the explicit formulation is typically favoured over implicit (i.e. strongly coupled) approaches for the fluid-structure interaction is also because of its computationally inexpensive nature. The detailed explanation of the methodology is, however, beyond the scope of this paper. For an in-depth discussion of the formulation, we refer readers to the paper of Spandan et al. (2017).
Numerical set-up for code validation consisting of a mixed convection flow over a rigid heated sphere. Here, U ∞ and θ ∞ denote the laminar free stream velocity and temperature, respectively. In addition, θ sph denotes sphere temperature, which is greater than θ ∞ . Also shown is a typical temperature field for Re sph = 100, Ri sph = 0.2 and Pr = 0.72, where red regions represent the normalised temperature value of one and yellow regions represent zero. Only a subset of the domain is shown. (b) Plot of Nu sph and C D versus increasing ratio of sphere diameter to local grid size, D/Δx. Percentages shown are relative to lower D/Δx values.

Code validation
The IBM code used in this study has been previously validated for various particle-flow configurations (e.g. de Tullio & Pascazio 2016; Spandan et al. 2018a;Chong et al. 2020). Given that the present study is a more complex flow system involving more parameters, here, we provide details of additional validation tests for a mixed convection problem. Specifically, our test set-up is a laminar flow over a rigid sphere which is held at a constant temperature (θ sph ) and is hotter than the free stream temperature (θ ∞ ), see figure 2(a). Here, the gravity vector opposes the free stream velocity and the resulting buoyancy flux is 'assisting' the flow (Kotouč, Bouchet & Dušek 2009;Musong & Feng 2014). The flow is periodic in y-and z-directions and an outflow, radiation boundary condition is prescribed at the outlet for velocities and temperature. The sphere is positioned in the middle of the domain.
For this test set-up, the governing parameters are the sphere Reynolds number (Re sph ≡ U ∞ D/ν), sphere Richardson number (Ri sph ≡ gβ(θ sph − θ ∞ )D/U 2 ∞ ) and Pr. The system responses are the sphere Nusselt number, Nu sph , and the drag coefficient,C D , defined as can be directly computed by numerical integration of the heat fluxes and stresses over the volume of the sphere. However, the numerical integration involves extending a probe normal to each discretised surface and performing additional MLS interpolation of velocities, pressure and temperature at the tip of the probe. Correspondingly, the number of calculations increases dramatically with increasingly finer resolutions (which are required for the following test cases), rendering the test simulations infeasible. The probe extension approach is also unnecessary since here we are dealing with a stationary, non-deforming mesh. Therefore, instead of (2.8a,b), we employ a more straightforward approach by directly integrating the immersed boundary thermal and hydrodynamic forcing terms over the N E Eulerian grid points associated with the Lagrangians, i.e. where (Breugem 2012;Musong & Feng 2014;Wang et al. 2019). Here, ΔV k is the forcing volume associated with each Eulerian point and is equal to the Eulerian cell volume. Hereinafter, we report Nu sph and C D values computed according to (2.9a,b).
First, we test the convergence for this set-up. For this test, the ratio of sphere diameter to local grid size (D/Δx) is varied from 16 to 80 and we fix Re sph = 60, Pr = 0.72, Ri sph = 4.0 and the domain size (L x × L y × L z = 8D × 4D × 4D). As discussed in § 2.1, the ratio D/Δx is a crucial parameter in IBM since a sufficiently small grid size is necessary to properly resolve the thermal and viscous boundary layers around the sphere (Wang et al. 2019). Figure 2(b) shows the trend of Nu sph and C D versus D/Δx. Also shown in the figure are the percentage of change relative to lower D/Δx values. From the figure, Nu sph and C D converges to within 0.3 % and 4.2 %, respectively, for D/Δx 48. In this particular test case C D appears to be more sensitive, however, as D/Δx is increased, the percentage change reduces to below 0.5 %. For our simulations for droplets in VC, D/Δx of 40 and higher are achieved in near-wall regions where the grid spacings are finer. Therefore, we conclude that the resolutions used in the VC flow with droplets are sufficiently resolved and a reasonable compromise in order to keep our simulations tractable.
Next, we validate our simulations with results from the literature. For these tests, we employed a larger domain where L x × L y × L z = 10D × 5D × 5D, in accordance with the domain sensitivity studies in Musong & Feng (2014). We vary Re sph (=60, 100), Ri sph (=0.2, 0.3, 4.0) and Pr (=0.72, 7.0). For this test, a more judicious numerical setting is warranted and, therefore, we used larger D/Δx values, where D/Δx ≈ 50 for Re sph = 60 and D/Δx ≈ 80 for Re sph = 100. The results are summarised in table 2.
From the table, our results for C D generally agree with results from the literature, especially for the Re sph = 100 cases. We note large differences for our C D values and the values of Wong, Lee & Chen (1986). Therefore, as an added assurance, we repeated the test case of the flow over a sphere at Re sph = 100 without the active temperature field (i.e. Ri sph = 0) and find that our results remain in good agreement to within 3.7 % when compared with simulations of Liska & Colonius (2017), Wang & Zhang (2011) and Kim, Kim & Choi (2001). For Nu sph , we find that our overall results are also in good agreement with the literature to within ≈4 % maximum. The consistency of our results indicate that our IBM achieves reasonable accuracy provided the grid is fine enough to resolve the boundary layers. However, since this approach inherently comes with a high computational cost, we choose to balance our IBM resolution strategy with the resolution necessary to resolve the background VC flow, as rationalised in § 2.1.

Description of the Lagrangian governing equations
Following the definition of IBM for deformable interfaces/fluid-structure interaction, the droplet interface is represented by a network of Lagrangian nodes evolved by the IP model (de Tullio & Pascazio 2016;Spandan et al. 2017). The equation of motion for each Lagrangian node, l, moving with velocity u l is (2.10) In (2.10), the terms are made dimensionless with the unit length L z and free-fall velocity U . The forces contributing to the right-hand side of (2.10) are the hydrodynamic loads F h , buoyancy F g and internal forces F i , where In (2.11a), V l is the volume of the node, but lacks a physical definition because the definition of the thickness of a liquid-liquid interface is not straightforward. To overcome this, following Spandan et al. (2017), we treat V l as a free parameter and fix V l = 1. Now, F h can be computed from direct integration of the viscous and pressure stresses from the flow, whereas F g is prescribed by varying Ra and Ra d (the definition for Ra d is described in § 2.5). The internal forces F i , on the other hand, come from modelling the droplet deformation characteristics using the IP model. The idea behind the IP model is to employ a spring network of nodes with tunable model parameters in order to represent the discretised droplet surface. On this point, a word of caution is warranted: because this is a model based on elastic structures, the IP model is an approximation of the actual interfacial dynamics arising from surface tension phenomena. The reason why the model is viable is because it is a phenomenological model, i.e. exact interfacial dynamics and the IP model both rely on the fundamental principle of minimum potential energy. However, the limitation of the model is that it cannot handle extreme deformations such as droplet breakup phenomenon because (i) the determination of the model parameters for large We is non-trivial, and (ii) the Lagrangian resolutions become terminally high.
A brief description of the IP model is as follows: F i represents the surface forces acting on the nodes of the discretised droplet surface. Under external hydrodynamic loads, the network of nodes deform and stores potential energy into the system. The potential energy is subsequently converted to surface forces by differentiating the potentials with respect to the displacements of each node. Details of the individual potentials of the IP model are outlined in § 2.3 of Spandan et al. (2017).

Model for thermally coupled droplets
For the lumped-capacitance model, two simplifying assumptions are made: (i) the droplets do not generate heat, and (ii) the internal temperature fields (and therefore interfacial temperature) of the droplets are uniform. Based on these assumptions, the interfacial droplet temperature is updated at every time step according to where θ d is the mean interfacial droplet temperature, V d is the volume of the droplet, S d is the droplet surface area and n is the outwardly directed unit normal. The droplet surface temperature is initialised as the mean surface temperature at its injected location, according to the initialisation step described earlier in § 2.1. After injection, the droplets rise and deform with respect to their original state (a sphere with diameter D), but do not significantly change in volume. Our model is therefore simpler than other numerical models with thermal coupling, for instance, studies that consider droplet growth at the boiling limit (e.g. Oresta et al. 2009;Lakkaraju et al. 2011) or models that rely on droplets with a constant geometry (e.g. Wang et al. 2017). The specific heat capacity ratio of the droplets to liquid phase (c p,d /c p,f ) is set equal to 2, so that the heat capacity ratios of the droplets to liquid phase (c d /c f ) are approximately equal to α. The assumption is based on the following: the total heat capacity of a phase is fixed by the specific heat, density and volume of the phase. Taking the ratio of heat capacities, we obtain c where c p is the specific heat capacity. Finally, with c p,d /c p,f = 2, which is equal toρ −1 in our set-up, c d /c f ≈ α. The heat capacity ratio therefore only becomes important at large α.

Derivation of the droplet Rayleigh number
In addition to the control parameters defined in (2.4a,b), we introduce the droplet Rayleigh number, Ra d , to quantify the droplet driving. It is defined as which is conveniently derived from scaling arguments of the governing equations outlined in § 2.3. We focus on the second term F g in (2.11b). Since (2.11b) represents the contribution from an isolated droplet and we are interested in defining a parameter for collective droplet effects, it would be reasonable to include the volume fraction parameter, α. Therefore, for 0 <ρ < 1, we define (2.14) which quantifies the relative dominance of droplet driving to thermal driving. It is important to keep in mind that in deriving (2.14), we did not consider other parameters (such as droplet size, heat capacity ratio and deformability) which would play a role towards the final flow dynamics (Serizawa, Kataoka & Michiyoshi 1975;Shen et al. 2006;Kitagawa & Murai 2013;Verschoof et al. 2016). In this study, we propose that (2.14) is informative when used as an engineering estimate to quantify similar flow problems. However, extensions to other flow configurations would require practical fine tuning based on systematic data, which are currently lacking. Other dimensionless parameters similar to Ra d /Ra have also been proposed for different flow configurations, but these require a priori knowledge of the dispersed phase dynamics and/or flow statistics. For example, Climent & Magnaudet (1999) proposed the Rayleigh number expression, Ra CM ≡ ρgαH 3 /(νU b ) (H is the height of the liquid layer and U b is the relative rise velocity of the bubble), to quantify bubble-induced convection. Based on the notion of pseudo-turbulence (Lance & Bataille 1991), which is defined as the fluctuating energy induced by the passage of bubbles under non-turbulent conditions, van Wijngaarden (1998) proposed the so-called bubblance parameter b ≡ (1/2)U 2 b α/u 2 0 (u 0 is the vertical velocity fluctuations of background turbulence). Since Ra d /Ra is a natural control parameter for VC with light droplets, we therefore use this ratio as input for our simulations. Note that Ra d is constant for a given α and therefore Ra d /Ra reduces with increasing Ra (this is equivalent to an increase in Froude number with increasing Ra). To make the simulations of the fluid-structure interaction tractable, we also run the simulations at g/200. The resulting Ra d /Ra is 5 × 10 −4 -5 × 10 −5 for α = 5 × 10 −3 and 2 × 10 −3 -2 × 10 −4 for α = 2 × 10 −2 .

Droplet influence on flow statistics and profiles
In this section, we analyse the results for 0 α 2 × 10 −2 , starting with a discussion of the droplets statistics.

Distribution of droplet aspect ratio versus bubble Reynolds number
From our simulations, the maximum droplet Reynolds number is Re d ≈ 220 and its time-averaged value is, Re d ≈ 100. As the droplets rise, they undergo deformation from the interfacial hydrodynamic loads. In figure 3, we characterise the deformation of the droplets in our simulations using the aspect ratio, Γ , of the horizontal to vertical axes (most often identical to the ratio of major to minor axes), which are determined by fitting two-dimensional Fourier descriptors (Duineveld 1995;Lunde & Perkins 1998) to the projected droplet outlines in the x y-and xz-plane. The joint probability density distribution in figure 3 shows that the droplets undergo moderate deformation between Γ ≈ 1 to Γ ≈ 1.3, agreeing with the relatively small We values. Values of Γ < 1 are  caused by small deformations in the droplet shapes after the reinjection step at the lower boundary, where the droplets are stirred by the cold fluid front. Visual inspections of the instantaneous shapes (insets of figure 3) show that the spherical droplet loses its fore-and-aft symmetry, with the front of the droplet becoming flatter than the back. Due to the relatively moderate Re d values, we do not observe droplet path instabilities throughout our simulations. In figure 3(b), we show that the droplets concentration profile is uniformly distributed for each respective α value and averaged over all Ra cases.

Profiles of mean vertical velocity and temperature
Now, we turn our focus to the flow statistics. To establish a baseline, we first analyse the influence of the droplets on the mean flow profiles of VC. Figure 4 shows the mean vertical velocity and temperature profiles plotted versus the horizontal component z (note that all length scales have been made dimensionless with L z ). Without droplets, the mean profiles are anti-symmetric about the channel centreline ( figure 4a,d). The cell centre is stably stratified (figure 6d) with d θ xy /dz| z=0.5 = 0 and u xy | z=0.5 = 0. Therefore, unlike the doubly periodic VC set-up (Ng et al. 2015(Ng et al. , 2017, there is no persistent mean shear in the bulk of the flow. For α > 0 and for both coupling cases, the mean vertical velocity profiles are asymmetric with a much stronger downward velocity magnitude near the cooler walls ( figure 4b,c). This symmetry-breaking phenomenon can be more clearly observed in figure 5, where the y-and time-averaged vertical velocity and temperature fields for Ra = 0.1 × 10 9 and α = 5 × 10 −3 are visualised. In figure 5(a), the downwards flowing (colder) fluid extends almost the entire vertical extent x as compared with the upwards flowing (hotter) fluid. The difference between the maxima and minima of u xy is largest for the smallest Ra, indicating that the droplet forcing is strongest. The mean temperature profiles (figure 4e, f ) also exhibit asymmetries. For α = 5 × 10 −3 and at the lowest Ra, the temperature profiles for both coupling cases are relatively constant and do not exhibit any undershoot, which is observed for α = 0 in figure 4(d) at z ≈ 0.04. However, at higher Ra, the profiles now bear some resemblance to the cases when α = 0, corroborating the notion that thermal driving increasingly dominates. Here, we note that although θ xy > θ ref in the bulk, the globally averaged temperature field θ is statistically stationary within 0.5 % for all cases. In the bulk region (0.2 z 0.8), we obtain d θ xy /dz| bulk ≈ 0. Based on these results, the influence of the light droplets is seemingly most pronounced at the vertical boundaries as compared to the bulk.  thermal boundary layers at the vertical walls. Since the blockage factor is higher for the α = 2 × 10 −2 cases, the magnitude of the mean horizontal velocities are larger at x 0.3 as compared to the α = 5 × 10 −3 cases. For the temperature profiles, we note an overall weakening of the stable stratification at higher α (figures 6e and 6f ), with the bulk mean temperatures θ yz → θ ref . The relatively uniform value of θ yz for the most part of x indicates strong mixing of the thermal field with increasing α.

Root-mean-square profiles of vertical velocity and temperature
The root mean square (r.m.s.) of the fluctuating quantities are plotted in figure 7 for all cases as a function of horizontal component z. Here, we define (·) rms ≡ ( u 2 xy ) 1/2 . When α = 0 (figure 7a,d), both u rms and θ rms exhibit near-wall peaks and are symmetrical about the channel centreline.
When α > 0, the bulk velocity fluctuations u bulk,rms > 0 as a direct result of droplet induced liquid fluctuations. Interestingly, u bulk,rms at lower Ra values are much larger than at higher Ra, which highlights the greater influence of droplet forcing on the flow at lower Ra. A rough estimate of the amplitude of the droplet-induced liquid agitations is as follows: for the lowest Ra, computing the ratios of max[u rms /U ] between α > 0 and α = 0 gives relative perturbation magnitudes of ≈ 3 and 6 times, for α = 5 × 10 −3 Increasing Ra and 2 × 10 −2 , respectively. The ratios are smaller at higher Ra because the background VC flow increasingly dominates the droplet-induced liquid agitations. The u rms profiles also exhibit slight asymmetry with values tending to be larger closer to the colder wall as compared to the hotter wall. This asymmetry is consistent with the notion of a more intermittent colder downwards flow caused by the disruption of the large-scale circulation by the droplet passage, as discussed in § 3.2. For θ rms , the magnitudes in the bulk for α > 0 (figure 7e, f ) tend to be lower than for the case when α = 0 (figure 7d), where θ rms,bulk ≈ 0.2. With thermal coupling, the θ rms profiles are typically slightly larger than without thermal coupling and counteracts the mechanical agitation by the droplets. This effect can be explained by the thermal exchange of the droplet and the surrounding liquid which induces local thermal fluctuations. Therefore, both the mechanical agitation at larger α and the thermal coupling of the droplets contribute to the bulk mixing of the thermal field.

Scaling of Nusselt and Reynolds numbers versus Rayleigh number
In figure 8, we present the scaling of the Nu and Re versus Ra. Here, we employ the wind-based Reynolds number, Re ≡ u xy,max L z /ν as a measure of the large-scale circulation.
When α = 0.0 (solid circles, figure 8), we find that Nu ∼ Ra 0.25±0.003 and Re ∼ Ra 0.50±0.002 which are in agreement with the Nu ∼ Ra 1/4 and Re ∼ Ra 1/2 analytical predictions for laminar boundary layer-dominated VC (Shishkina 2016). For VC with droplets (triangles), the Nu values appear to be larger at higher Ra values; however, a big variation about their mean persists across the range of Ra simulated. Given this uncertainty, it is therefore unclear whether a power law exists in the present parameter space and we dispense with any attempts to fit effective power laws. Further judicious studies at larger separations of Ra values would be prudent and could provide more information. On the other hand, the Re trends are clearly less steep than the Ra 1/2 scaling with Re values that are much larger at lower Ra and decrease in magnitude with increasing Ra. When comparing the coupling cases, the effective scaling for Nu and Re is largely unaffected. However, by including thermal coupling, the temperature field is distributed more efficiently, and so the magnitude of the heat transport is increased. Albeit small, this increase is an interesting result because our small Biot number assumption implies a weaker influence of thermal coupling in our flow. As a direct comparison for Nu, the ratio Nu/Nu α=0 is shown in the inset of figure 8(a) and the values range from 0.95 to 1.1. Some caution is warranted here when interpreting the ratios. Because of the rather large variations of NuRa −1/4 as shown in the figure, we cannot conclusively claim that there exists a decrease in Nu at low Ra. However, we can link the variations of the ratios to the different manner in which the droplets locally influence the wall heat fluxes and wall shear stresses. The local influences are quantified and discussed in § 4 and in § 5. Now, we focus on the Re trends. For α > 0, the Re values tend to be larger than for the α = 0 case and this is consistent with the response of the VC flow due to the passage of the droplets across the top and bottom boundary layers. As the droplets cross the horizontal boundary layers, the large-scale circulation of the background VC flow is continuously disrupted, triggering horizontal intrusions of warmer fluid at the top wall and cooler fluid at the bottom wall (peaks in mean horizontal velocities in figures 6b and 6c), similar to the intrusions observed in transient VC in a square cavity (Patterson & Imberger 1980;Armfield & Patterson 1991). For α = 5 × 10 −3 , at the higher Ra values, the Re values tend to approach the Re values for α = 0. This incipient trend suggests that the droplet driving is no longer dominant at this part of the parameter space as compared to the α = 2 × 10 −2 case.

Droplet influence on local Nusselt number
In this section, we link the Nu versus Ra variations discussed in § 3.5 to the changes in the local Nusselt number evaluated at the hot and cold walls. We define the local Nusselt number as Nu loc ≡ f w,loc L z /(Δκ) = [∂ θ y (x)/∂z]| w /( /L z ), which is the local dimensionless temperature gradient evaluated at z = 0 and L z . The trends are shown in figure 9 as function of x.
From figure 9, Nu loc are larger in the upstream of the vertical boundary layers, that is x 1.2 for figure 9(a-c) and x 1.2 for figure 9(d-f ). Here, the larger values of Nu loc simply reflect the thinner thermal boundary layers developing from the corners of the domain. For α = 0, Nu loc monotonically decreases as the boundary layer develops and is consistent across the Ra range. However, the trends vary considerably for α > 0. For example, relative to the α = 0 cases, (i) Nu loc,h becomes lower for x 1.2, and (ii) for α = 2 × 10 −2 , both Nu loc,h and Nu loc,c are roughly constant for 0.6 x 1.8. Since these changes directly reflect the thermal boundary layer thicknesses, we can conclude that the droplets not only influence the bulk statistics as shown in § 3, but would also influence the local thermal boundary layers.
To Mech. + therm. coupling behaviour can also be observed for α = 2 × 10 −2 , although the corresponding wall area with decreased Nu loc,h is smaller for the mechanically coupled case (see figure 10c and the inset plot). The decreased Nu loc,h for the α = 5 × 10 −3 case overwhelms the increased Nu loc,c for x 1.2, with the lowest Ra cases being most strongly influenced, as previously shown in figure 8. In contrast, Nu loc,c is significantly increased for α = 2 × 10 −2 and x 1.2 by roughly a factor of 1.5 times (figure 10d,h). Based on the much stronger droplet driving for α = 2 × 10 −2 , Nu| α=2×10 −2 is increased by approximately 5 % for the lowest Ra relative to Nu| α=5×10 −3 . For the different distributions of Nu loc for α > 0 in figures 9 and 10, we note that the profiles represent local quantities and are non-monotonic at both walls with increasing Ra. The spatial distributions therefore cannot be trivially determined a priori. What can be discerned from the current results is that the droplets influence the bulk flow (as seen in the mean and r.m.s. statistics in figures 4 to 7), the near-wall flow and the large-scale circulation of VC. Different mechanisms in these regions compete and the prevailing mechanism(s) would presumably determine the heat transport of the set-up.

Droplet influence on local skin-friction coefficient
Unlike Rayleigh-Bénard convection, the thermal boundary layers in VC are sheared by a mean wind with a constant direction that is predetermined by the boundaries (Ng et al. 2015). Therefore, to quantify the influence of the droplets on wind shearing, we plot the local skin-friction coefficient C f ,loc versus x in figure 11. Here, C f ,loc ≡ 2τ w (x)/U 2 , where . Similar to figure 9, but now for the local skin-friction coefficient C f ,loc . Colour legends are the same as figure 4.
τ w (x) ≡ μ∂ u y (x)/∂z| w is the wall shear stress. Similar to the idea of figure 10, the relative changes in the local skin-friction coefficients are plotted in figure 12. For α = 0, C f ,loc is largest at wall heights that are close to the upstream of the developing boundary layer. However, when α > 0, C f ,loc is roughly constant for the most part of x at low Ra. Two points can be made from the distributions of C f ,loc . First, the roughly uniform distribution of C f ,loc at low Ra for α > 0 imply that the droplet driving dominates the mean wind of VC and, on a mean sense, homogenises the viscous boundary layer particularly at the hot wall. Second, the distributions of C f ,loc are not symmetric at the hot and cold walls (for example, max [C f ,loc,c ] > max [C f ,loc,h ]) as compared to the α = 0 case (figures 10a and 10d). One possible explanation of this asymmetry can be made by observing the rising direction of the droplets: at the cold wall, the droplets oppose the downwards flow whereas at the hot wall, the droplets aid the upwards flow. Coupled with the asymmetry of the mean horizontal velocity profiles in figure 6, the resulting viscous boundary layer becomes thinner at the cold wall, and a larger C f ,loc results. However, this conjecture may not hold at higher Ra cases because the viscous boundary layers eventually become much thinner and closer to the walls. As a result, at sufficiently high Ra, the influence of droplets presumably diminishes with increasing distance from the edge of the viscous boundary layers, eventually yielding to the dynamics of thermal driving.
When compared with C f ,loc,0 (figure 12), we find larger values of C f ,loc in concomitant regions with larger values of Nu loc in figure 10. Interestingly, whilst Nu loc is relatively insensitive to Ra (see figure 10), the wall-height distributions of C f ,loc exhibit a strong non-monotonic behaviour which depends on Ra, α and whether the cold or hot wall  Mech. + therm. coupling is considered. Therefore, it appears that C f ,loc is more sensitive to the droplets induced agitation as compared to Nu loc .
6. Light droplets versus bubbles -a comparison to experiments by Gvozdić et al. (2018) In this section, we discuss several aspects of the physical parameters in our simulations, which distinguish our findings from the laboratory results of bubbly VC by Gvozdić et al. (2018).
A crucial difference between our investigation and the experiments is thatρ = 0.5 in our simulations (corresponding to light droplets) whereasρ ≈ 10 −3 in their experiments (corresponding to air bubbles in water). Clearly, the large differences in the density ratios play a role and this is reflected in our simulations. For example, the mean temperature in the bulk region of our simulations have approximately zero gradient (figure 4), whereas the mean temperature in the bulk region of the experiments have a finite gradient (see figure  9(a) of their paper), indicating a much stronger mixing of the thermal field by the bubbles as compared to light droplets. Furthermore, the values of Nu for VC with light droplets is within 10 % of the Nu values without droplets (figure 8), whereas in the laboratory experiments of Gvozdić et al. (2018), Nu can be larger by up to 20 times with bubbles than without and remain Ra-independent for their investigated parameter range. Therefore, we conclude that the background VC flow remains relatively dominant even with influence of light droplets, and this is reflected in the non-monotonic distributions of the local heat transport and skin-friction coefficients shown in figures 10 and 12.
The strength of the bubble-induced agitation versus droplet-inducted agitation can also be quantified a posteriori using the bubblance parameter b ≡ U 2 b α/u 2 0 (6.1) (cf. Lance & Bataille 1991;van Wijngaarden 1998;Rensen, Luther & Lohse 2005;Alméras et al. 2017), which defines the ratio of energy produced by a bubble swarm, i.e. U 2 b α, and the energy of the background turbulence without bubbles, i.e. u 2 0 . Note that a prefactor of one is chosen for (6.1), which is different to previous definitions which employ a prefactor of 1/2 (based on the added mass coefficient, cf. Rensen et al. (2005)); however, the present discussions are still valid. We define U b as the mean bubble or droplet rise velocity and u 0 as the maximum of the mean vertical velocity of the single phase flow at the half-height of the domain. Next, we estimate b for our DNS and for the experiments by Gvozdić et al. (2018).
Since we have the full information from our DNS, the calculation of b is straightforward. For the laboratory experiments, u 0 was not recorded and so, invoking dynamic similarity, we estimate the values using our DNS results at matched Ra. Here, U b is assumed to be 0. It is useful for applications such as in chemical mixing, to have an estimate of the parameter space for b or Ra d /Ra where the driving by background turbulence eventually dominates bubble driving. For the laboratory experiments with bubbly VC, Gvozdić et al. (2018) estimated this parameter space by defining a critical Rayleigh number, Ra c , as follows: first, an effective power law trend of Nu ∼ Ra 0.33 is obtained from the single phase experiments. Then, observing that the Nu trends are insensitive to Ra for 5 × 10 −3 α 5 × 10 −2 (cf. figure 12 of Gvozdić et al. (2018)), the Nu ∼ Ra 0.33 and constant Nu trends are extrapolated to higher Ra values. The intersection of these curves is defined as Ra c , where 7 × 10 10 Ra c (α) 2 × 10 12 for the α values investigated. The range of Ra c values is marked in figure 13(c).
We can now directly extrapolate the trends of b and Ra d /Ra to the Ra c values. From least square fits, the effective power laws are b ∼ Ra −1 and Ra d /Ra ∼ Ra −1 . Therefore, the extrapolated values are b ∼ Ra −1 c and Ra d /Ra ∼ Ra −1 c , visually marked by the blue patches in figure 13. For illustration purposes, only the α = 5 × 10 −3 and 5 × 10 −2 are drawn and an allowance of Ra c ± 10 % was employed to compute the extrapolation. The corresponding values are (b, Ra d /Ra)| α=5×10 −3 ≈ (0.2, 60) and (b, Ra d /Ra)| α=5×10 −2 ≈ (0.06, 18). These values suggest that the VC flow will dominate bubble-induced liquid agitation at b O(10 −1 ) and Ra d /Ra O(100). We note that our dataset for α = 2 × 10 −2 coincide with this regime for b| α=5×10 −2 (lower horizontal blue line in figure 13a), however, since the boundary layer dynamics are still dominant for our configuration, it suggests thatρ is an additionally important parameter when characterising bubbly turbulence. Interestingly, for bubbles rising in grid-generated turbulence (or incident turbulence), Alméras et al. (2017) determined a slightly larger value for b (≈0.7), where bubble-induced agitation appears to dominate. The mechanism was related to an increase in development length of the secondary bubble wake, which significantly enhances liquid velocity fluctuations. Indeed, the values of b from our DNS are smaller which is consistent with the notion that the background flow remains dominant for our parameter space considered.

Conclusions and outlook
In this study, we simulated the VC flow with dispersed light droplets between Ra = 1.3 × 10 8 and 1.3 × 10 9 and Pr value of 7. The liquid phase is simulated using DNS whereas the dispersed phase is simulated using an IBM with the interaction potential method for deformable interfaces. Our approach extends the IBM of Spandan et al. (2017) and Spandan et al. (2018a), where now the dispersed phase is fully coupled both mechanically and thermally to the flow. In addition, two datasets are simulated with and without thermal coupling to investigate its influence on the heat transport. Although Nu is slightly larger when the droplets are thermally coupled, we found that the VC flow with light droplets exhibits a non-monotonic change in heat transport with increasing Ra and largely retains the laminar-like VC scaling. We reason that a significant enhancement of heat transport depends crucially on a sufficiently strong droplet driving, which we show can be characterised by the relative strength of Ra d to Ra and the bubblance parameter, b.
When light droplets are introduced, the mean velocity and temperature profiles are highly skewed with the lowest Ra being most sensitive (figures 4-7). However, this sensitivity is masked by the Nu versus Ra trend, where we observe a non-monotonic behaviour with increasing Ra (figure 8a). This suggests the presence of competing mechanisms in the flow that contribute to the net heat transport. In contrast, the decreasing Re versus Ra trends are commensurate with the higher sensitivity at lower Ra, i.e. mechanical stirring is strongest at lowest Ra and higher α (figure 8b).
Based on analyses of the near-wall regions, we found that regions with higher values of local heat fluxes, Nu loc , correspond to concomitant regions with higher values of skin-friction coefficient, C f ,loc , which is consistent with the notion that the local wind has influence over the local heat transport (figures 9 and 11). In turn, the strength of the local wind is related to whether the direction of the rising droplets aids or opposes the flow ( figure 12). However, the trends of Nu loc and C f ,loc remain spatially non-monotonic and is sensitive to α for the simulation parameters considered in this study.
The Nu versus Ra trend in figure 8 is different from recent experimental results by Gvozdić et al. (2018) for bubbly flow. Whilst Nu exhibits some Ra-dependency for our simulations with light droplets, Gvozdić et al. (2018) reported that Nu is largely insensitive to Ra for various volume fractions of droplets. The key distinction between our DNS and the experiments by Gvozdić et al. (2018) becomes readily apparent when we quantify the bubblance parameter b and the droplet driving parameter Ra d /Ra (cf. § 6). Both b and Ra d /Ra have a large separation in scales between the laboratory experiments and our DNS. More specifically, at b O(10 −1 ) and Ra d /Ra O(100), we anticipate that the dynamics of the dispersed phase-induced liquid agitations become overwhelmed by the dynamics of the background VC flow. For light droplets, both b and Ra d /Ra are significantly lower. Therefore the local heat fluxes and skin friction coefficients exhibit non-monotonic behaviour, which reflects the dominance of the background VC flow.
Our results collectively indicate a non-monotonic heat transport behaviour for light droplets. Locally, the near-wall trends of heat fluxes and wall-shear stresses suggest the presence of competing mechanisms that, in concert, govern heat transport. One question that arises naturally here is: Can monotonicity be eventually obtained by increasing b and Ra d /Ra for fixed Ra? The answer to this question may provide some clues about disentangling the competing heat transport mechanisms in multiphase VC and warrants judicious numerical studies at larger parameter spaces.