2.1. 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 Reference Verzicco and Orlandi1996). 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:

1. (i) the resolution at the corner flow regions;

2. (ii) the resolution at the bulk flow; and

3. (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. Reference Shishkina, Stevens, Grossmann and Lohse2010). As a check, a coarser simulation with 20 % fewer grid points results in $\textit{Nu}$ values that are within 0.5 %, indicating good convergence for our resolution. For point (ii), we determined that $\max [{\rm \Delta} x^+_i \equiv {\rm \Delta} x_i/\delta _\nu ]\approx 2.4$ (details in table 1), where ${\rm \Delta} x_i$ are the grid spacings in each $i$th-direction and $\delta _\nu \equiv \nu /u_\tau$ is the viscous length scale based on the local shear velocity scale $u_\tau \equiv [\nu (\partial \langle u \rangle _{y}(x)/\partial z)|_w]^{1/2}$. Here, $\langle \cdot \rangle _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/(\text {max}[{\rm \Delta} x_i]) \gtrapprox 28$, comparable to other immersed boundary studies in turbulent flow with finite-size particles (Wang, Vanella & Balaras Reference Wang, Vanella and Balaras2019). Other numerical strategies are certainly possible, such as employing uniform grid spacings (Lu, Fernández & Tryggvason Reference Lu, Fernández and Tryggvason2005) or by eliminating walls in the simulations (Uhlmann & Chouippe Reference Uhlmann and Chouippe2017), 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,

(2.7)\begin{equation} {\textit{Re}}_d \equiv U_dD/\nu, \end{equation}

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 Reference Johnson and Patel1999). Given that our maximum droplet rise Reynolds number, $\max [{\textit {Re}}_d]\approx 220$, for an equivalent sphere Reynolds number, the dimensionless boundary layer thickness at its stagnation point is $\delta _{sp}/D\approx 1.13{\textit {Re}}_d^{-1/2} \approx 0.08$ (Schlichting & Gersten Reference Schlichting and Gersten2000). 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 ${\textit {Re}}_d$ value and largest grid spacing in our set-up. Our IBM resolution improves at lower ${\textit {Re}}_d$ (i.e. for thicker droplet boundary layers) and for finer near-wall grid spacings.

Table 1. Summary of simulation parameters. The corresponding number of grid points are $(n_x,n_y,n_z) = (960,96,384)$ for ${\textit {Ra}} \leqslant 0.4\times 10^{9}$ and $(n_x,n_y,n_z) = (1200,120,480)$ for ${\textit {Ra}} \geqslant 0.7\times 10^{9}$. Here, $T_{s}/(L_z/U_\Delta )$ and $T_s/\langle t_{d}\rangle$ are the total simulation sampling interval in terms of the free-fall velocity and droplet rise times, respectively.

Our IBM employs the fast MLS algorithm (Spandan et al. Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017, Reference Spandan, Lohse, de Tullio and Verzicco2018a). Two volume fractions are simulated: $\alpha = 5\times 10^{-3}$ and $2\times 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 $\alpha = 5\times 10^{-3}$ and 48 droplets for $\alpha = 2\times 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. Reference Gvozdić, Alméras, Mathai, Zhu, van Gils, Verzicco, Huisman, Sun and Lohse2018) where the ratio of bubble diameter to channel width ${\sim } O(10^{-2})$. Sizes of dispersed phase have been shown to play a role in influencing momentum and heat transport (Shen, Ceccio & Perlin Reference Shen, Ceccio and Perlin2006; Kitagawa & Murai Reference Kitagawa and Murai2013; Verschoof et al. Reference Verschoof, van Der Veen, Sun and Lohse2016); 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 ${\textit {Re}}_d \approx O(100)$ of our simulations is $L_w/D \approx 1$ (cf. Clift, Grace & Weber (Reference Clift, Grace and Weber2005), Chapter 5). Given that $L_y/D \approx 3.1 > 1$, the spanwise domain size is sufficiently large and we assume that the droplets do not interact with their own wake.

At the start of each simulation, droplets are spatially initialised as spheres in a $4\times 3$ array (in the $xz$-plane at $y=0.5L_y$) for $\alpha = 5\times 10^{-3}$ and a $12\times 4$ array for $\alpha = 2\times 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, $\Delta /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/\langle t_{d}\rangle$ where $\langle t_{d}\rangle$ 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 ${\textit {Ra}}$ (defined by $T_{s}/(L_z/U_\Delta )$, where $U_\Delta \equiv (g\beta {\rm \Delta} 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 $\alpha = 2\times 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. Reference Spandan, Verzicco and Lohse2018b). Nevertheless, we still employed the collision detection algorithm of Spandan et al. (Reference Spandan, Lohse, de Tullio and Verzicco2018a) in our numerics and the elastic potential collision model by Spandan et al. (Reference Spandan, Verzicco and Lohse2018b). 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, $\theta _{{init}}^k$, is equated to the averaged temperature of the immersed fluid, according to $\theta _{{init}}^k = N_f^{-1} \sum _{k=1}^{N_f} \theta ^k$, where $\theta ^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; $\theta _{{init}}^k$ is subsequently forced using IBM to the Eulerian grid. A small initial droplet vertical velocity ${\sim } O(10^{-2})U_\Delta$ 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 ${\textit {Bi}} \equiv 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 Reference Wang, Sierakowski and Prosperetti2017). Our reasoning for the small droplet Biot number value is as follows: by substituting the length scale $l_d \equiv D/6$ (ratio of sphere volume to sphere surface area) and the heat transfer coefficient $h \equiv {\textit {Nu}}_d k_f/D$, we obtain ${\textit {Bi}} = {\textit {Nu}}_d k_f/(6k_d)$. Here, ${\textit {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 ${\textit {Nu}}_d \sim O(1)$ (cf. Oresta et al. Reference Oresta, Verzicco, Lohse and Prosperetti2009). Finally, we assume $k_d > k_f$, and so ${\textit {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 Reference Duineveld1995; Jenny, Dušek & Bouchet Reference Jenny, Dušek and Bouchet2004). 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, ${\textit {We}} \equiv \rho _{{ref}} U_\Delta ^2 D/\sigma$, which yields the ratio of inertia to capillary forces, where $\sigma$ is the surface tension. In our simulation strategy, $\sigma$ is not prescribed explicitly. Rather, an additional tuning step is performed to obtain a set of IP model parameters such that ${\textit {We}} \approx 3\times 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. (Reference Spandan, Verzicco and Lohse2018b). 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 Reference Maffettone and Minale1998) and a droplet in cross-flow (e.g. Loth Reference Loth2008; Schwarz, Kempe & Fröhlich Reference Schwarz, Kempe and Fröhlich2016). Finally, in accordance with the tuning criteria described in Spandan et al. (Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017), $\sigma$ is reverse-engineered by matching the droplet deformation dynamics to the reference results in Maffettone & Minale (Reference Maffettone and Minale1998) and Schwarz et al. (Reference Schwarz, Kempe and Fröhlich2016). 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 ${\textit {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 $\hat {\rho }$, $D/L_z$ or $\alpha$. After extensive precursor simulations and checks, we decided to simulate droplets at half the density of the fluid, i.e. $\hat {\rho }\equiv \rho _d/\rho _{{ref}} = 0.5$, which is within the numerical stability limit for explicit IBM time integration schemes (Schwarz, Kempe & Fröhlich Reference Schwarz, Kempe and Fröhlich2015). The value of $\hat {\rho }$ 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. (Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017).

2.2. Code validation

The IBM code used in this study has been previously validated for various particle-flow configurations (e.g. de Tullio & Pascazio Reference de Tullio and Pascazio2016; Spandan et al. Reference Spandan, Lohse, de Tullio and Verzicco2018a; Chong et al. Reference Chong, Li, Ng, Verzicco and Lohse2020). 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 ($\theta _{{sph}}$) and is hotter than the free stream temperature ($\theta _\infty$), 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 Reference Kotouč, Bouchet and Dušek2009; Musong & Feng Reference Musong and Feng2014). 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.

Figure 2. ($a$) Numerical set-up for code validation consisting of a mixed convection flow over a rigid heated sphere. Here, $U_\infty$ and $\theta _\infty$ denote the laminar free stream velocity and temperature, respectively. In addition, $\theta _{{sph}}$ denotes sphere temperature, which is greater than $\theta _\infty$. Also shown is a typical temperature field for ${\textit {Re}}_{{sph}}=100$, ${\textit {Ri}}_{{sph}}=0.2$ and $\textit {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 ${\textit {Nu}}_{{sph}}$ and $C_D$ versus increasing ratio of sphere diameter to local grid size, $D/{\rm \Delta} x$. Percentages shown are relative to lower $D/{\rm \Delta} x$ values.

For this test set-up, the governing parameters are the sphere Reynolds number (${\textit {Re}}_{{sph}} \equiv U_\infty D/\nu$), sphere Richardson number (${\textit {Ri}}_{{sph}} \equiv g\beta (\theta _{{sph}}-\theta _\infty )D/U_\infty ^2$) and $\textit {Pr}$. The system responses are the sphere Nusselt number, $\widetilde {{\textit {Nu}}}_{{sph}}$, and the drag coefficient, $\tilde {C}_D$, defined as

(2.8a,b)\begin{equation} \widetilde{{\textit{Nu}}}_{{sph}} \equiv \widetilde{F_\theta}/\left[{\rm \pi} (\theta_{{sph}}-\theta_\infty) D\right] \quad\text{and}\quad \tilde{C}_D \equiv \widetilde{F_x}/\left[(1/2)\rho_\infty U_\infty^2{\rm \pi}(D/2)^2\right], \end{equation}

where $\widetilde {F_\theta } = -\oint _{\partial V}\boldsymbol {\nabla }\theta \boldsymbol {\cdot }\boldsymbol {n}\,\textrm {d} A$ and $\widetilde {F_x} = \oint _{\partial V}\delta _{i1}(\boldsymbol {\tau }\boldsymbol {\cdot }\boldsymbol {n})\,\textrm {d} A$. Equation (2.8a,b) 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.

(2.9a,b)\begin{equation} {\textit{Nu}}_{{sph}} \equiv F_\theta/\left[{\rm \pi} (\theta_{{sph}}-\theta_\infty) D\right] \quad\text{and}\quad C_D \equiv F_x/\left[(1/2)\rho_\infty U_\infty^2{\rm \pi}(D/2)^2\right], \end{equation}

where $F_\theta \equiv -\kappa ^{-1}\sum _{k=1}^{N_{E,{tot}}}f^k_\theta \,{\rm \Delta} V^k$ and $F_x \equiv -\rho _\infty \sum _{k=1}^{N_{E,{tot}}}f^k_x\,{\rm \Delta} V^k$ (Breugem Reference Breugem2012; Musong & Feng Reference Musong and Feng2014; Wang et al. Reference Wang, Vanella and Balaras2019). Here, ${\rm \Delta} V^k$ is the forcing volume associated with each Eulerian point and is equal to the Eulerian cell volume. Hereinafter, we report ${\textit {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/{\rm \Delta} x$) is varied from 16 to 80 and we fix ${\textit {Re}}_{{sph}}=60$, $\textit {Pr}=0.72$, ${\textit {Ri}}_{{sph}}=4.0$ and the domain size ($L_x\times L_y\times L_z = 8D \times 4D \times 4D$). As discussed in § 2.1, the ratio $D/{\rm \Delta} 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. Reference Wang, Vanella and Balaras2019). Figure 2(b) shows the trend of ${\textit {Nu}}_{{sph}}$ and $C_D$ versus $D/{\rm \Delta} x$. Also shown in the figure are the percentage of change relative to lower $D/{\rm \Delta} x$ values. From the figure, ${\textit {Nu}}_{{sph}}$ and $C_D$ converges to within $0.3\,\%$ and $4.2\,\%$, respectively, for $D/{\rm \Delta} x \geqslant 48$. In this particular test case $C_D$ appears to be more sensitive, however, as $D/{\rm \Delta} x$ is increased, the percentage change reduces to below 0.5 %. For our simulations for droplets in VC, $D/{\rm \Delta} 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\times L_y\times L_z = 10D \times 5D \times 5D$, in accordance with the domain sensitivity studies in Musong & Feng (Reference Musong and Feng2014). We vary ${\textit {Re}}_{{sph}}$ (${=}60, 100$), ${\textit {Ri}}_{{sph}}$ (${=}0.2, 0.3, 4.0$) and $\textit {Pr}$ (${=}0.72, 7.0$). For this test, a more judicious numerical setting is warranted and, therefore, we used larger $D/{\rm \Delta} x$ values, where $D/{\rm \Delta} x \approx 50$ for ${\textit {Re}}_{{sph}} = 60$ and $D/{\rm \Delta} x \approx 80$ for ${\textit {Re}}_{{sph}} = 100$. The results are summarised in table 2.

Table 2. Lift coefficients $C_D$ and sphere Nusselt numbers ${\textit {Nu}}_{{sph}}$ at various ${\textit {Re}}_{{sph}}$, ${\textit {Ri}}_{{sph}}$ and $\textit {Pr}$. Also provided are results from MF14 – Musong & Feng (Reference Musong and Feng2014); KB09 – Kotouč et al. (Reference Kotouč, Bouchet and Dušek2009); WL86 – Wong et al. (Reference Wong, Lee and Chen1986); LC17 – Liska & Colonius (Reference Liska and Colonius2017); WZ11 – Wang & Zhang (Reference Wang and Zhang2011); and KK01 – Kim et al. (Reference Kim, Kim and Choi2001).

$^{a}$Interpolated from data.

From the table, our results for $C_D$ generally agree with results from the literature, especially for the ${\textit {Re}}_{{sph}} = 100$ cases. We note large differences for our $C_D$ values and the values of Wong, Lee & Chen (Reference Wong, Lee and Chen1986). Therefore, as an added assurance, we repeated the test case of the flow over a sphere at ${\textit {Re}}_{{sph}} = 100$ without the active temperature field (i.e. ${\textit {Ri}}_{{sph}}=0$) and find that our results remain in good agreement to within $3.7\,\%$ when compared with simulations of Liska & Colonius (Reference Liska and Colonius2017), Wang & Zhang (Reference Wang and Zhang2011) and Kim, Kim & Choi (Reference Kim, Kim and Choi2001). For ${\textit {Nu}}_{{sph}}$, we find that our overall results are also in good agreement with the literature to within ${\approx }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.

2.3. 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 Reference de Tullio and Pascazio2016; Spandan et al. Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017). The equation of motion for each Lagrangian node, $l$, moving with velocity $\boldsymbol {u}_l$ is

(2.10)\begin{equation} \dfrac{\textrm{d} \boldsymbol{u}_l}{\textrm{d} t} = \boldsymbol{F}^{h} + \boldsymbol{F}^{g} + \boldsymbol{F}^{i}. \end{equation}

In (2.10), the terms are made dimensionless with the unit length $L_z$ and free-fall velocity $U_\Delta$. The forces contributing to the right-hand side of (2.10) are the hydrodynamic loads $\boldsymbol {F}^{h}$, buoyancy $\boldsymbol {F}^{g}$ and internal forces $\boldsymbol {F}^{i}$, where

(2.11a,b)\begin{equation} {\boldsymbol{F}}^{h} = \dfrac{L_z}{\hat{\rho}V_lU_\Delta^2} \int_{S} \tau \boldsymbol{\cdot} {\boldsymbol{n}}\,\textrm{d} S \quad \text{and} \quad {\boldsymbol{F}}^{g} \equiv \left(1-\dfrac{1}{\hat{\rho}}\right)\dfrac{L_z}{U_\Delta^2} {\boldsymbol{g}}. \end{equation}

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. (Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017), we treat $V_l$ as a free parameter and fix $V_l=1$.

Now, $\boldsymbol {F}^{h}$ can be computed from direct integration of the viscous and pressure stresses from the flow, whereas $\boldsymbol {F}^{g}$ is prescribed by varying ${\textit {Ra}}$ and ${\textit {Ra}}_d$ (the definition for ${\textit {Ra}}_d$ is described in § 2.5). The internal forces $\boldsymbol {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 ${\textit {We}}$ is non-trivial, and (ii) the Lagrangian resolutions become terminally high.

A brief description of the IP model is as follows: $\boldsymbol {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. (Reference Spandan, Meschini, Ostilla-Mónico, Lohse, Querzoli, de Tullio and Verzicco2017).

2.4. 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

(2.12)\begin{equation} \dfrac{\textrm{d} \theta_d}{\textrm{d} t} ={-}\dfrac{\kappa}{V_d} \oint_{S_d} \boldsymbol{\nabla} \theta \boldsymbol{\cdot} \boldsymbol{n}\,\textrm{d} S_d, \end{equation}

where $\theta _d$ is the mean interfacial droplet temperature, $V_d$ is the volume of the droplet, $S_d$ is the droplet surface area and $\boldsymbol {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. Reference Oresta, Verzicco, Lohse and Prosperetti2009; Lakkaraju et al. Reference Lakkaraju, Schmidt, Oresta, Toschi, Verzicco, Lohse and Prosperetti2011) or models that rely on droplets with a constant geometry (e.g. Wang et al. Reference Wang, Sierakowski and Prosperetti2017).

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 $\alpha$. 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_d/c_f = (c_{p,d}/c_{p,f})\hat {\rho }\alpha /(1-\alpha ) \approx (c_{p,d}/c_{p,f})\hat {\rho }\alpha$ for small $\alpha$ values, where $c_p$ is the specific heat capacity. Finally, with $c_{p,d}/c_{p,f} = 2$, which is equal to $\hat {\rho }^{-1}$ in our set-up, $c_d/c_f \approx \alpha$. The heat capacity ratio therefore only becomes important at large $\alpha$.

2.5. Derivation of the droplet Rayleigh number

In addition to the control parameters defined in (2.4a,b), we introduce the droplet Rayleigh number, ${\textit {Ra}}_d$, to quantify the droplet driving. It is defined as

(2.13)\begin{equation} {\textit{Ra}}_d \equiv \dfrac{\alpha gL_z^3}{\hat{\rho}\nu\kappa}, \end{equation}

which is conveniently derived from scaling arguments of the governing equations outlined in § 2.3.

We focus on the second term $\boldsymbol {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, $\alpha$. Therefore, for $0 < \hat {\rho } < 1$, we define

(2.14)\begin{equation} F^{g}_\alpha \sim \dfrac{\alpha gL_z}{\hat{\rho}U_\Delta^2} =: \dfrac{{\textit{Ra}}_d}{{\textit{Ra}}}, \end{equation}

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 Reference Serizawa, Kataoka and Michiyoshi1975; Shen et al. Reference Shen, Ceccio and Perlin2006; Kitagawa & Murai Reference Kitagawa and Murai2013; Verschoof et al. Reference Verschoof, van Der Veen, Sun and Lohse2016). 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 ${\textit {Ra}}_d/{\textit {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 (Reference Climent and Magnaudet1999) proposed the Rayleigh number expression, ${\textit {Ra}}_{CM} \equiv \rho g \alpha H^3/(\nu 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 Reference Lance and Bataille1991), which is defined as the fluctuating energy induced by the passage of bubbles under non-turbulent conditions, van Wijngaarden (Reference van Wijngaarden1998) proposed the so-called bubblance parameter $b \equiv (1/2)U_b^2\alpha /u_0^{2}$ ($u_0$ is the vertical velocity fluctuations of background turbulence). Since ${\textit {Ra}}_d/{\textit {Ra}}$ is a natural control parameter for VC with light droplets, we therefore use this ratio as input for our simulations. Note that ${\textit {Ra}}_d$ is constant for a given $\alpha$ and therefore ${\textit {Ra}}_d/{\textit {Ra}}$ reduces with increasing ${\textit {Ra}}$ (this is equivalent to an increase in Froude number with increasing ${\textit {Ra}}$). To make the simulations of the fluid–structure interaction tractable, we also run the simulations at $g/200$. The resulting ${\textit {Ra}}_d/{\textit {Ra}}$ is $5\times 10^{-4}$–$5\times 10^{-5}$ for $\alpha = 5\times 10^{-3}$ and $2\times 10^{-3}$–$2\times 10^{-4}$ for $\alpha = 2\times 10^{-2}$.

3.1. Distribution of droplet aspect ratio versus bubble Reynolds number

From our simulations, the maximum droplet Reynolds number is ${\textit {Re}}_d \approx 220$ and its time-averaged value is, $\langle {\textit {Re}}_d\rangle \approx 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, $\varGamma$, 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 Reference Duineveld1995; Lunde & Perkins Reference Lunde and Perkins1998) to the projected droplet outlines in the $xy$- and $xz$-plane. The joint probability density distribution in figure 3 shows that the droplets undergo moderate deformation between $\varGamma \approx 1$ to $\varGamma \approx 1.3$, agreeing with the relatively small ${\textit {We}}$ values. Values of $\varGamma < 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 ${\textit {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 $\alpha$ value and averaged over all ${\textit {Ra}}$ cases.

Figure 3. ($a$) Joint probability density function of droplet deformations characterised by $\varGamma$ versus the droplet rise Reynolds number, ${\textit {Re}}_d$, averaged over all cases. Outermost contour level is 0.03 and the contours are spaced 0.06 apart (reproduced in the colourbar for emphasis). Inset plots of panel ($a$) show representative two-dimensional droplet shapes for two ${\textit {Re}}_d$ values. Also shown are the directions of the droplet rise velocity $U_d$ and gravity $g$. ($b$) Normalised droplet concentration profiles for $\alpha = 5\times 10^{-3}$ (blue) and $2\times 10^{-2}$ (red) as a function of horizontal component $z$, averaged across all ${\textit {Ra}}$. The associated colour-shaded regions show ${\pm } 1\sigma$ variation about the averaged concentration value at the corresponding $z$ location.

3.2. 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 $\textrm {d} {\langle \theta \rangle _{xy}}/\textrm {d} z|_{z=0.5} = 0$ and ${\langle u \rangle _{xy}}|_{z=0.5}=0$. Therefore, unlike the doubly periodic VC set-up (Ng et al. Reference Ng, Ooi, Lohse and Chung2015, Reference Ng, Ooi, Lohse and Chung2017), there is no persistent mean shear in the bulk of the flow. For $\alpha > 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 ${\textit {Ra}} = 0.1\times 10^9$ and $\alpha = 5\times 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 ${\langle u \rangle _{xy}}$ is largest for the smallest ${\textit {Ra}}$, indicating that the droplet forcing is strongest.

Figure 4. Mean profiles as a function of horizontal component $z$: ($a$–$c$) vertical velocity ${\langle u\rangle _{xy}}/U_\Delta$; ($d$–$f$) temperature $({\langle \theta \rangle _{xy}}-\theta _{{ref}})/\Delta$. ($a,d$) $\alpha =0$; ($b,e$) $\alpha =5\times 10^{-3}$; and ($c,f$) $\alpha =2\times 10^{-2}$. For $\alpha =0$, only the left half of the profiles are shown since the profiles are antisymmetric about the vertical centreline. Dashed grey curves represent mechanical coupling only. Solid red curves represent mechanical and thermal coupling. Darker curves represent higher ${\textit {Ra}}$.

Figure 5. Visualisations of (a) vertical velocity, and (b) temperature field for ${\textit {Ra}} = 0.1\times 10^9$ and $\alpha = 5\times 10^{-3}$. The streamlines depict the cellular-like large-scale circulation and the distortions from the passage of droplets. In panel (b), the upstream and downstream regions are indicated at the hotter (red) and cooler (blue) walls, respectively.

Figure 6. Same as figure 4, but now for the mean profiles as a function of the vertical component $x$: (a–c) horizontal velocity ${\langle w\rangle _{yz}}/U_\Delta$; (d–f) temperature $({\langle \theta \rangle _{yz}}-\theta _{{ref}})/\Delta$. For $\alpha = 0$, only the velocity profiles between $0\leqslant x \leqslant 1.2$ are shown since the profiles are antisymmetric about the horizontal centreline. Colour legends are the same as figure 4.

The mean temperature profiles (figure 4e,f) also exhibit asymmetries. For $\alpha = 5\times 10^{-3}$ and at the lowest ${\textit {Ra}}$, the temperature profiles for both coupling cases are relatively constant and do not exhibit any undershoot, which is observed for $\alpha = 0$ in figure 4(d) at $z\approx 0.04$. However, at higher ${\textit {Ra}}$, the profiles now bear some resemblance to the cases when $\alpha = 0$, corroborating the notion that thermal driving increasingly dominates. Here, we note that although ${\langle \theta \rangle _{xy}} > \theta _{{ref}}$ in the bulk, the globally averaged temperature field ${\langle \theta \rangle }$ is statistically stationary within 0.5 % for all cases. In the bulk region ($0.2 \leqslant z \leqslant 0.8$), we obtain $\textrm {d} {\langle \theta \rangle _{xy}}/\textrm {d} z|_{{bulk}}\approx 0$. Based on these results, the influence of the light droplets is seemingly most pronounced at the vertical boundaries as compared to the bulk.

3.3. Profiles of mean horizontal velocity and temperature

Figure 6 shows the mean horizontal velocity and temperature profiles plotted versus the vertical component $x$. When $\alpha = 0$, the velocity profiles are antisymmetric (figure 6a) and the temperature profiles are constant for all ${\textit {Ra}}$ (figure 6d). When $\alpha > 0$, the antisymmetries are destroyed: for $\alpha = 5 \times 10^{-3}$, the horizontal velocities are larger at the top wall (figure 6b), whereas for $\alpha = 2\times 10^{-2}$, the horizontal velocities are larger at the bottom wall (figure 6c). The source for the asymmetry can be traced to the passage of droplets entering the bottom or leaving the top of the domain: at the lower boundary, the droplets which have near-zero velocity block the horizontal flow causing the fluid to accelerate around the droplets. At the upper boundary, the droplets exit the domain at terminal velocity, and the entrained fluid impinges on the upper wall. Both mechanisms trigger intermittent intrusions of hotter and colder fluid at the upstream corners of the thermal boundary layers at the vertical walls. Since the blockage factor is higher for the $\alpha = 2\times 10^{-2}$ cases, the magnitude of the mean horizontal velocities are larger at $x \lesssim 0.3$ as compared to the $\alpha = 5\times 10^{-3}$ cases.

For the temperature profiles, we note an overall weakening of the stable stratification at higher $\alpha$ (figures 6e and 6f), with the bulk mean temperatures ${\langle \theta \rangle _{yz}} \rightarrow \theta _{{ref}}$. The relatively uniform value of ${\langle \theta \rangle _{yz}}$ for the most part of $x$ indicates strong mixing of the thermal field with increasing $\alpha$.

3.4. 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 $(\cdot )'_{rms}\equiv (\langle u'\rangle _{xy}^2)^{1/2}$. When $\alpha = 0$ (figure 7a,d), both $u'_{rms}$ and $\theta '_{rms}$ exhibit near-wall peaks and are symmetrical about the channel centreline.

Figure 7. Same as figure 4, but now for the r.m.s. statistics as a function of the horizontal component $z$: (a–c) $10u'_{rms}/U_\Delta$; (d–f) $\theta '_{rms}/\Delta$. For $\alpha = 0$, only the left half of the profiles are shown since they are antisymmetric about the vertical centreline. Colour legends are the same as figure 4.

When $\alpha > 0$, the bulk velocity fluctuations $u'_{{bulk},rms} > 0$ as a direct result of droplet induced liquid fluctuations. Interestingly, $u'_{{bulk},rms}$ at lower ${\textit {Ra}}$ values are much larger than at higher ${\textit {Ra}}$, which highlights the greater influence of droplet forcing on the flow at lower ${\textit {Ra}}$. A rough estimate of the amplitude of the droplet-induced liquid agitations is as follows: for the lowest ${\textit {Ra}}$, computing the ratios of max[$u'_{rms}/U_\Delta$] between $\alpha > 0$ and $\alpha = 0$ gives relative perturbation magnitudes of $\approx 3$ and $6$ times, for $\alpha = 5\times 10^{-3}$ and $2\times 10^{-2}$, respectively. The ratios are smaller at higher ${\textit {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 $\theta '_{{rms}}$, the magnitudes in the bulk for $\alpha > 0$ (figure 7e,f) tend to be lower than for the case when $\alpha = 0$ (figure 7d), where $\theta '_{{rms},{bulk}} \approx 0.2$. With thermal coupling, the $\theta '_{{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 $\alpha$ and the thermal coupling of the droplets contribute to the bulk mixing of the thermal field.

3.5. Scaling of Nusselt and Reynolds numbers versus Rayleigh number

In figure 8, we present the scaling of the ${\textit {Nu}}$ and ${\textit {Re}}$ versus ${\textit {Ra}}$. Here, we employ the wind-based Reynolds number, ${\textit {Re}} \equiv \langle u\rangle _{xy,{max}}L_z/\nu$ as a measure of the large-scale circulation.

Figure 8. (a) Compensated ${\textit {Nu}}$ versus ${\textit {Ra}}$, and (b) compensated ${\textit {Re}}$ versus ${\textit {Ra}}$. Inset of panel (a): ratio of ${\textit {Nu}}|_{\alpha > 0}$ to ${\textit {Nu}}|_{\alpha = 0}$. Solid black symbols are DNS data for $\alpha = 0$. Upwards-pointing blue triangles are for $\alpha = 5\times 10^{-3}$ cases and downwards-pointing red triangles are for $\alpha = 2\times 10^{-2}$. Open triangles denote mechanical coupling only and filled triangles denote both mechanical and thermal coupling. The ${\textit {Nu}}$ versus ${\textit {Ra}}$, and ${\textit {Re}}$ versus ${\textit {Ra}}$ scalings for $\alpha = 0$, are consistent with analytical predictions for VC driven by laminar boundary layers, i.e. ${\textit {Nu}}\sim {\textit {Ra}}^{1/4}$ and ${\textit {Re}}\sim {\textit {Ra}}^{1/2}$ at constant $\textit {Pr}$ (Shishkina Reference Shishkina2016).

When $\alpha = 0.0$ (solid circles, figure 8), we find that ${\textit {Nu}} \sim {\textit {Ra}}^{0.25\pm 0.003}$ and ${\textit {Re}} \sim {\textit {Ra}}^{0.50\pm 0.002}$ which are in agreement with the ${\textit {Nu}} \sim {\textit {Ra}}^{1/4}$ and ${\textit {Re}} \sim {\textit {Ra}}^{1/2}$ analytical predictions for laminar boundary layer-dominated VC (Shishkina Reference Shishkina2016). For VC with droplets (triangles), the ${\textit {Nu}}$ values appear to be larger at higher ${\textit {Ra}}$ values; however, a big variation about their mean persists across the range of ${\textit {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 ${\textit {Ra}}$ values would be prudent and could provide more information. On the other hand, the ${\textit {Re}}$ trends are clearly less steep than the ${\textit {Ra}}^{1/2}$ scaling with ${\textit {Re}}$ values that are much larger at lower ${\textit {Ra}}$ and decrease in magnitude with increasing ${\textit {Ra}}$. When comparing the coupling cases, the effective scaling for ${\textit {Nu}}$ and ${\textit {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 ${\textit {Nu}}$, the ratio ${\textit {Nu}}/{\textit {Nu}}_{\alpha =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 ${\textit {Nu}}{\textit {Ra}}^{-1/4}$ as shown in the figure, we cannot conclusively claim that there exists a decrease in ${\textit {Nu}}$ at low ${\textit {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 ${\textit {Re}}$ trends. For $\alpha > 0$, the ${\textit {Re}}$ values tend to be larger than for the $\alpha = 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 Reference Patterson and Imberger1980; Armfield & Patterson Reference Armfield and Patterson1991). For $\alpha = 5\times 10^{-3}$, at the higher ${\textit {Ra}}$ values, the ${\textit {Re}}$ values tend to approach the ${\textit {Re}}$ values for $\alpha = 0$. This incipient trend suggests that the droplet driving is no longer dominant at this part of the parameter space as compared to the $\alpha = 2\times 10^{-2}$ case.