Experimental investigation of heat transport in homogeneous bubbly flow

We present results on the global and local characterisation of heat transport in homogeneous bubbly flow. Experimental measurements were performed with and without the injection of $\sim 2.5$ mm diameter bubbles (corresponding to $Re_b \approx 600$) in a rectangular water column heated from one side and cooled from the other. The gas volume fraction $\alpha$ was varied in the range $0\% - 5\%$, and the Rayleigh number $Ra_H$ in the range $4.0 \times 10^9 - 1.2 \times 10^{11}$. We find that the global heat transfer is enhanced up to 20 times due to bubble injection. Interestingly, for bubbly flow, for our lowest concentration $\alpha = 0.5\% $ onwards, the Nusselt number $\overline{Nu}$ is nearly independent of $Ra_H$, and depends solely on the gas volume fraction~$\alpha$. We observe the scaling $\overline{Nu} \propto \alpha^{0.45}$, which is suggestive of a diffusive transport mechanism. Through local temperature measurements, we show that the bubbles induce a huge increase in the strength of liquid temperature fluctuations, e.g. by a factor of 200 for $\alpha = 0.9\%$. Further, we compare the power spectra of the temperature fluctuations for the single- and two-phase cases. In the single-phase cases, most of the spectral power of the temperature fluctuations is concentrated in the large-scale rolls. However, with the injection of bubbles, we observe intense fluctuations over a wide range of scales, extending up to very high frequencies. Thus, while in the single-phase flow the thermal boundary layers control the heat transport, once the bubbles are injected, the bubble-induced liquid agitation governs the process from a very small bubble concentration onwards.

For example, in vertical natural convection, the usage of fins or riblets is proven to increase the heat flux (Shakerin et al. 1988). In Rayleigh-Bénard convection Lohse & Xia 2010) heat transfer enhancement can be achieved by tuning the boundary conditions to aid the formation of large scale rolls or coherent structures (Chong et al. 2015), or by introducing wall roughness elements (e.g. Roche et al. (2001); Tisserand et al. (2011); Xie & Xia (2017); Zhu et al. (2017b)). While these methods have been widely used to optimize convective transport, they pose limits on the maximum achievable heat flux in moderate to high Rayleigh number flows.
An alternative is the injection of bubbles in the flow, which indeed enhances the heat transport considerably (Deckwer 1980). In general, the bubbles can be injected in a quiescent liquid phase ("pseudo-turbulent" flow (Risso & Ellingsen 2002;Riboux et al. 2010;Roghair et al. 2011;Mercado Martínez et al. 2010)) or in an already turbulent liquid phase (turbulent bubbly flow (Rensen et al. 2005;Van Den Berg et al. 2006;van Gils et al. 2013;Spandan et al. 2016;Prakash et al. 2016)). It is known that the motion of injected bubbles induces mixing of warm and cold parcels of the liquid phase, which in industrial applications where heat transport is coupled with the bubbly flow can lead to a 100 times greater heat transfer coefficient when compared to the single-phase case (Deckwer 1980). Therefore there is a practical benefit from a better understanding of the heat transport in bubbly flows as this enables better design and optimization of the industrial processes. As a result, the effect of bubbles on heat transfer has been subject of several experimental and numerical studies in the past.
The bubbles can be injected in a system with natural or forced convection. Early studies which focused on forced convective heat transfer in bubbly flows (Sekoguchi et al. 1980;Sato et al. 1981a,b) showed that the bubbles modify the temperature profile and that higher void fractions close to the heated wall lead to an enhanced heat transfer. In a more recent numerical study on forced convective heat transfer in turbulent bubbly flow in vertical channels, Dabiri & Tryggvason (2015) showed that both, nearly spherical and deformable bubbles, improve the heat transfer rate. They found that a 3% volume fraction of bubbles increases the Nusselt number by 60%.
Studies on natural convection in bubbly flow have been performed mostly by introducing micro-bubbles (Kitagawa et al. 2008(Kitagawa et al. , 2009) and sub-millimeter-bubbles (Kitagawa & Murai 2013) close to the heated wall. Among these, the micro-bubbles (mean bubble diameter d bub = 0.04 mm) showed higher heat transfer enhancement as compared to sub-millimeter-bubbles (d bub = 0.5 mm) (Kitagawa & Murai 2013). The authors stated that this occurred because the micro-bubbles form large bubble swarms which rise close to the wall and enhance mixing in the direction of the temperature gradient, while submillimeter bubbles, owing to their weak wake and low bubble number density, resulted in limited mixing. In contrast, in case of bubbles with diameter of a few millimeters, the wake of individual bubbles and vortex shedding behind the bubbles play a significant role in the heat transfer enhancement (Kitagawa & Murai 2013). Tokuhiro & Lykoudis (1994) experimentally studied the effects of 2 − 4 mm diameter inhomogeneously injected nitrogen-bubbles on laminar and turbulent natural convection heat transfer from a vertical heated plate in mercury. They reported a twofold increase in the heat transfer coefficient as compared to the case without bubbles. Similarly, Deen & Kuipers (2013) showed in their numerical study that a few high Reynolds number bubbles rising in quiescent liquid could increase the local heat transfer between the liquid and a hot wall.
In systems with natural convection bubbles can be introduced through boiling as well. Some of the studies on this subject have been performed for Rayleigh-Bénard convection. The Rayleigh-Bénard system consists of a flow confined between two horizontal parallel plates, where the bottom plate is heated and the top one is cooled Lohse & Xia 2010). In case of boiling it has been found that bubbles strongly affect velocity and temperature fields depending on the Jakob number which is defined as a ratio of latent heat to sensible heat (Oresta et al. 2009;Zhong et al. 2009;Schmidt et al. 2011). Numerical studies performed by Lakkaraju et al. (2013) found that without taking into consideration the bubble nucleation and bubble detachment, depending on the number of the bubbles and superheat the heat transfer can be enhanced up to around 6 times.
In an attempt to control the bubble nucleation process, Narezo Guzman et al. (2016b,a) performed an experimental study where they varied the geometry of the nucleation sites of the bubbles and found that the heat transfer could be enhanced up to 50%.
To summarize, previous studies (performed in systems with natural convection) mostly focused either on heat transfer in inhomogeneous bubbly flow, where bubbles of different sizes were introduced close to the hot wall, and where the thermal stability of the used setups remained unclear, or the studies were performed in a well-defined Rayleigh-Bénard system, but in which the bubbles mainly consisted of vapor and where bubble volume fraction and the bubble size varied due to evaporation and condensation.
In this study on heat transfer in bubbly flow we choose a different approach. Firstly, we use a rectangular water column heated from one side and cooled from the other (see figure 1) which resembles the classical vertical natural convection system. Secondly, at the bottom of the setup we homogeneously inject millimetric-bubbles, so that in the bubbly case pseudo-turbulent flow is present, the dynamics of which is adequately characterized and broadly studied in the past. We characterise the global heat transfer in both the single-and two-phase flow cases with the goal of understanding how the imposed temperature difference (characterised by the the Rayleigh number Ra H ) influences the heat flux (characterised by the Nusselt number N u) for various gas volume fractions α, in order to try to better understand the mechanism of heat transport enhancement. The characterization of the global heat transfer is based on the calculation of the dimensionless temperature difference, the Rayleigh number: (1.1) and the dimensionless heat transfer rate, the Nusselt number: Here, Q is the measured power supplied to the heaters, T h and T c are the mean temperatures of the hot and cold walls, respectively, L is the length of the setup, A is the surface of the sidewall, β is the thermal expansion coefficient, g the gravitational acceleration, κ the thermal diffusivity, and χ the thermal conductivity of water. In this study we choose the height to be the characteristic length scale for Ra H since the boundary layer regime is present (see Section 3), where the velocity is predominately in the vertical direction. Note that in this study the Prandtl number is nearly constant, P r ≡ ν κ = 6.5 ± 0.3. Furthermore, we perform temperature profile measurements with and without bubble injection, by traversing a small thermistor along the length of the setup at the mid-height. In this way we obtain information on the statistics of the temperature fluctuations and on the power spectrum of the temperature fluctuations in a homogeneous bubbly flow.

Experimental setup
The experiments were performed in a rectangular bubble column (600×230×60 mm 3 ), shown in figure 1. Air bubbles of about 2.5 mm diameter were injected into quiescent demineralized water using 180 capillaries (inner diameter 0.21 mm) uniformly distributed over the bottom of the column. The gas volume fraction was varied from 0.5% to 5% by controlling the inlet gas flow rate via a digital mass flow controller (Bronkhorst F-111AC-50K-AAD-22-V). The two main sidewalls of the setup (600 × 230 mm 2 ) were made of 1 cm thick glass and two (heated resp. cooled) sidewalls (600 × 60 mm 2 ) of 1.3 cm thick brass.
As mentioned previously, our setup resembles one for vertical natural convection since one brass sidewall is heated and the other is cooled in order to generate a horizontal temperature gradient in the setup. More precisely, heating was provided by placing three etched-foil heaters on the outer side of the hot wall. Heaters were connected in parallel to a digitally controlled power supply (Keysight N8741A), providing altogether up to 300 W. The other brass wall was cooled by a water circulating bath (Polytemp PD15R-30). The temperature of the heated and the cooled walls was monitored by three thermistors which were glued on different heights of the cold and warm walls, namely at 125 mm, 315 mm, and 505 mm. Temperature regulation of both the cold and warm walls was achieved by PID (Proportional-Integral-Derivative) control so that the mean temperature of the walls was maintained constant over time. In order to limit the heat losses, the setup was wrapped in several layers of insulating blanket and foam. Moreover, an aluminium plate with heaters attached to it was placed on the outer side of the hot wall and maintained at the same temperature as the hot wall to act as a temperature shield.
Heat losses were estimated to be not more than 7% by calculating convective heat transport rate from all outer surfaces of the setup with the assumption that these surfaces are at maximum 25 • C. On the other hand, we measured the power needed to maintain the temperature of the bulk constant (T bulk = 25 • C) over 4 hours. Power supplied to the heaters in that way is not more that 3% of the total power needed when running the actual experiments in which the bulk temperature is also 25 degrees. We therefore expect the actual heat losses to be in the range of 3% to 7% which enables us to study precisely the heat transport driven by a horizontal temperature gradient in a bubbly flow.

Single phase vertical convection
The single-phase heat transport will be used as a reference for the heat transport enhancement by bubble injection. In a single-phase system a variety of flow regimes can be observed depending on the height H, the length L, and the Rayleigh number Ra H of the system (Bejan 2004). For the parameter range studied here (H/L = 2.4 and Ra H = 4.0 × 10 9 − 1.2 × 10 11 ) we expect the system to be in a boundary layer dominant regime. In this regime the boundary layers which distinctly form along the heated and cooled sidewalls control the heat transfer, while the bulk of the fluid is relatively stagnant. Furthermore, previous studies on single-phase vertical convection describe the dependence of Nusselt number on the Rayleigh number in the power law form N u ∼ Ra β at fixed P r, with exponent β ranging between 1/4 and 1/3 (see Ng et al. (2015Ng et al. ( , 2017 and references within). We find that for the range of Ra H studied here the effective scaling exponent is β ≈ 0.33 ± 0.02, which lies within the expected range as can be seen on the compensated plot in figure 2.
For the single phase flow we benchmark the heat transfer against direct numerical simulations (DNS). For this purpose an in-house second order finite difference code (van der Poel et al. 2015;Zhu et al. 2017a) was used to solve the three-dimensional Boussinesq equations for the single phase vertical convection. The code has been extensively validated and used for Rayleigh-Bénard flow (van der Poel et al. 2015;Zhu et al. 2017a), in which the only difference is the direction of the buoyancy force. The computational box has the same size as has been employed in the experiments. The no-slip boundary conditions are adopted for the velocity at all solid boundaries. At the top and bottom walls, the heat-insulating conditions are employed, and at the left and right plates, constant temperatures are prescribed. The resolution of the simulations is fine enough to guarantee that the results are grid-independent. Figure 2 shows good agreement between numerical and experimental results, within the range of uncertainty. As the experimental setup was not completely insulated at the top and due to unavoidable heat losses to the outside, the experimentally obtained N u are slightly higher.

Instrumentation for the gas phase characterization
The homogeneity of the bubble swarm was verified by measuring the gas volume fraction α at different locations within the flow by means of a single optical fiber probe. In the absence of a temperature gradient, we observe that α is nearly uniform along the length L and height H of the setup (see figure 3 (a)). We also note that no large-  scale clustering is present, since the cumulative distribution function F (∆t) of the time between consecutive bubbles ∆t is Poissonian (see figure 3 (b) and Risso & Ellingsen (2002) for more details). In presence of the heating, the homogeneity of the swarm was comparable to that without heating, with not more than 2% variation. The bubble swarm thus remained homogeneous even with heating, indicating that the bubbly flow was not destabilised by the temperature gradient. We further characterised the gas phase by measuring the bubble diameter and the bubble rising velocity with an in-house dual optical fiber probe (Alméras et al. 2017). Bubble diameter d eq and bubble rise velocity  (2014), we find up to 6% variation which is acceptable since it can be attributed to slightly different bubble injection section and water quality.

Instrumentation for the heat flux and temperature measurements
In the present study, we performed global and local characterisations of the heat transfer. In order to obtain N u and Ra H , we measured the hot and cold wall temperatures, and the heat input to the system Q. To this end, resistances of the thermistors placed on the hot and cold walls were read out every 4.2 seconds using a digital multimeter (Keysight 34970A). The temperature is then converted from the resistances based on the calibrations of the individual thermistors. The total heater power input was measured as where I i and V i are the current and voltage across each heater, respectively. The experimental measurements were performed after steady state was achieved in which the mean wall temperatures fluctuated less that ±0.01 K (resp. ±0.1 K) for lowest Ra H and ±0.4 K (resp. ±0.5 K) for highest Ra H , for single-phase (resp. two-phase) case. Time averaging of the instantaneous power supplied to each heater Q i was then performed over a total time period of 6 hours for single-phase cases, and 3 hours for two-phase cases.
For a better understanding of the heat transfer, we performed local measurements of the liquid temperature fluctuations: T = T − T , where T -temperature fluctuations, T -measured instantaneous temperature and T -time averaged temperature at the measurement point. For each operating condition, temperature fluctuations were recorded for at least 180 min (resp. 360 min) in two-phase (resp. single-phase) case once the flow was stable, which yields satisfactory statistical convergence. We used a NTC miniature thermistor manufactured by TE connectivity (Measurement Specialties G22K7MCD419) with a tip diameter of 0.38 mm, and a response time of 30 ms in water. The thermistor was connected as one arm of a Wheatstone bridge so that very small variations of the thermistors' resistance could be measured. For noise reduction we used a Lock-In amplifier (SR830). The function of the Lock-In amplifier is to firstly supply voltage to the bridge and then filter out noise at frequencies which are different from that of the signal range. This ensured milli-Kelvin resolution of the temperature fluctuations. In order to measure the temperature with this resolution reliably, the thermistors were calibrated in a circulating bath with 5 mK stability. A typical obtained signals in single -phase and two-phase are presented in the figure 5 a) and b), respectively. It is interesting to note that the temperature fluctuations are up to 200 times stronger in the two-phase case and this will be discussed in Section 4.
The local temperature measurement technique used here is well established for singlephase flow (Belmonte et al. 1994); however, until now it has never been used for temperature fluctuations measurements in bubbly flows. Since the presence of bubbles may perturb the local measurements by interacting directly with the probe (Rensen et al. 2005;Mercado Martínez et al. 2010), it is necessary to validate the technique for two-phase flow. For this purpose we made an in-house probe in which an optical fiber and the thermistor were positioned ∼ 1 mm apart in the horizontal plane and with the thermistor placed ∼ 1 mm below the fiber tip.  bubble-thermistor contact (t int = d eq /V bub = 7.2 ms for α = 0.9%) were filtered out due to the longer response time of the thermistor (t r ∼ 30 ms). Nevertheless, the thermistor response time is sufficiently short to measure the temperature fluctuations between two bubble passages. From table 1 we see that t 2b is sufficiently long when compared to t r even for higher gas volume fractions. Furthermore, the estimated time of the bubblethermistor contact remains much shorter than t r for α up to 5% (see table 1). The present measurement technique is thus suitable for measurements of the temperature fluctuations in the liquid phase not only for α = 0.9%, but also for higher gas volume fractions.

Global heat transport enhancement
Let us now discuss the heat transport in the presence of a homogeneous bubble swarm for gas volume fraction α ranging from 0.5% to 5%, and the Rayleigh number ranging from 4.0 × 10 9 to utmost 3.6 × 10 10 (see Figure 6). For the whole range of α and Ra H , adding bubbles considerably increases the heat transport, since the Nusselt number is about an order of magnitude higher as compared to single-phase flow (see figure 6 (a)). In order to better quantify the heat transport enhancement due to bubble injection, we show in figure 6 (b) the ratio of the Nusselt number in the bubbly flow N u b to the Nusselt number in the single-phase N u 0 as a function of Ra H for different α. We find that heat transfer is enhanced up to 20 times due to bubble injection, and that the enhancement increases with increasing α and decreasing Ra H . Note that the decreasing trend of N u b /N u 0 with Ra H occurs because the single-phase heat flux N u 0 increases with Ra H . Figure 6 (a) also shows that for a fixed α 0.5% N u remains nearly constant with increasing Ra H . This indicates that the boundary layers developing along the walls are not limiting the heat transport anymore in the two-phase case. Together with the observation that the heating does not induce a gradient in the gas volume fraction profile (see Figure 7), this further implies that the temperature behaves as a passive scalar in bubbly flow. In order to understand the mechanism of the heat transport in bubbly flow we compare the findings of our study to those of Alméras et al. (2015), who showed that the transport of a passive scalar at high Péclet number (P e = V bub d eq /D ≈ 10 6 , with V bub as the bubble velocity and d eq as the bubble diameter and D as molecular diffusivity) by a homogeneous bubble swarm is a diffusive process. In the case of a diffusive process, the turbulent heat flux can be modeled as u i T = −D ii ∇T , introducing the effective diffusivity D ii . If we take the heat transport to be a diffusive process in our study, the Nusselt number can be interpreted as the ratio between the effective diffusivity induced by the bubble swarm D ii and the thermal diffusivity χ. Thus, from our measurements we can estimate the effective diffusivity induced by a bubble swarm for a gas volume fraction ranging from 0.5% to 5%. Note that since the temperature gradient is imposed in the horizontal direction, we assume that the measured effective diffusivity is mainly in the horizontal direction. Figure 8 (a) shows the Nusselt number N u averaged over the full range of Rayleigh number for a constant gas volume fraction as a function of α. We clearly see that the averaged Nusselt number evolves as α 0.45±0.025 . Even if we subtract the single-phase Nusselt number (which might be thought of as the contribution of natural convection to the total heat transfer) from the one in bubbly flow the scaling remains unchanged (see figure 8 (b)). This trend is in a good agreement with the model of effective diffusivity proposed by Alméras et al. (2015). In fact, the authors showed that at low gas volume fraction, the diffusion coefficient can be written as: D ii ∝ u Λ, where u is the standard deviation of the velocity fluctuations, and Λ is the integral Lagrangian length scale (Λ d bub /C d0 , C d0 = 4d bub g/3V 0 2 , here C d0 is the drag coefficient and V 0 is the rise velocity of a single bubble (Riboux et al. 2010)). In the expression for the diffusion coefficient only the u depends on the gas volume fraction and this dependence is given by u ∼ V 0 α 0.4 (Risso & Ellingsen 2002). In the present study, we expect to have a similar liquid agitation since bubble rising velocity and diameter are comparable. This yields the same evolution of the effective diffusivity D ii with the gas volume fraction α namely, D ii ∝ α 0.45 , extending the model proposed by Alméras et al. (2015) to lower Péclet number (P e ≈ 5000 in this study). We also must stress here that the influence of the Péclet number on the effective diffusivity can be significant. In fact, numerical simulation performed by Loisy (2016), at α = 2.4% and Re b = 30 show that the effective diffusivity normalised by the molecular/thermal diffusivity varies linearly with the Péclet number (for P e ranging from 10 3 to 10 6 ). Consequently, since the Péclet number varies by three decades between the present study and the one of Alméras et al. (2015), no quantitative comparison of the effective diffusivity can be performed. Therefore, further studies on the effect of the Péclet number on the effective diffusivity at high Reynolds number should be performed.

Local characterisation of the heat transport
In order to gain further insight into the heat transport enhancement, we performed local liquid temperature measurements by traversing the thermistor along the length of the setup at mid-height for three Rayleigh numbers (Ra H = 5.2 × 10 9 , Ra H = 1.6 × 10 10 , and Ra H = 2.2 × 10 10 ), and for α = 0% and α = 0.9%. Figure 9 (a) shows the normalised temperature profiles. For the single-phase cases, as expected from the range of Ra H and the H/L in our study, a flat temperature profile along the length is observed in the bulk at mid-height (Elder 1965;Markatos & Pericleous 1984;Bejan 1984;Belmonte et al. 1994;Ng et al. 2015;Shishkina & Horn 2016;Ng et al. 2017). After normalisation, the single-phase temperature profiles for all three Rayleigh numbers overlap. However, due to present heat losses to the outside at the top of the setup since the setup is open on the top, the temperature profiles do not collapse at T −Tc T h −T c = 0.5 but at T −Tc T h −T c = 0.4. This has to be taken into account when comparing numerical data with the data obtained experimentally. Therefore, after shifting the numerically obtained temperature profiles good agreement is found between the two. The spatial temperature gradient in the singlephase case is located in the thermal boundary layer whose thickness δ t is estimated from the numerical data to be O(1 mm). In figure 10 we show the boundary layer thickness in the single phase case as a function of Ra H . Note that our flow configuration is different from the one present in classical Rayleigh-Bénard setup. Here the thermal boundary layer is defined as wall distance to the intercept of T = T h + dT /dx| w x and T = T h − ∆T /2 and the kinetic boundary layer is given as an intercept of u = du/dx| w x and u = u max (see e.g. Ng et al. (2015) for more details). The thickness of the thermal boundary layer based on the experimentally obtained Nusselt number is comparable to the one obtained numerically (see Figure 10).
As seen in figure 9 (a), the mean temperature profiles in the case of bubbly flow is completely different from that of single phase flow. In the bulk, two known mixing mechanisms contribute to the distortion of the flat temperature profile that was observed in the single-phase case: (i) capture and transport by the bubble wakes (Bouche et al. 2013), and (ii) dispersion by the bubble-induced turbulence which is the dominant one as shown by Alméras et al. (2015). Near the heated and cooled walls, we visually observed  (a)) and symbols (experimental data) present single-phase, blue symbols present two-phase with α = 0.9%, for various Rayleigh numbers: RaH = 5.2 × 10 9 (downwards triangles), RaH = 1.6 × 10 10 (squares), and RaH = 2.2 × 10 10 (upward triangles)).
the bubbles bouncing along the walls. This presumably disturbs the thermal boundary layers; however, this cannot be measured due to insufficient experimental resolution. Figures 9 (b) and 5 (a) and 5 (b) show that the temperature fluctuations induced by bubbles are even two orders of magnitudes higher than in the single-phase case (T = T − T , where T is the temperature fluctuations, T is the measured instantaneous temperature and T is the time averaged temperature at the measurement point). The normalised standard deviations of the fluctuations in both single-phase and two-phase are higher closer to the cold and hot walls than in the center of the setup. This is possibly due to the temperature probe seeing more hot and cold plumes closer to the heated and cold walls, a well known phenomenon from Rayleigh-Bénard flow ). Here slight asymmetry of the temperature profiles and the profiles of normalised standard deviation close to the walls must be attributed to the difference in the nature of heating and cooling of the sidewalls (see asymmetry also in figure 9 (a)). Figures 5 (a) and (b) also demonstrate that the time scales of fluctuations in single-phase and two-phase are different, along with much more intense temperature fluctuations for the bubbly flow.
To explore this in better detail, we now present the power spectrum of the temperature fluctuations ("thermal power spectrum") for both cases. Figure 11 (a) shows this power spectrum of temperature fluctuations at mid-height in the centre of the setup for the single-and two-phase cases. In the single-phase case the measured temperature fluctuations are limited to frequencies lower than 10 −1 Hz. At around 10 −2 Hz we observe a peak, beyond which there is a very steep decrease of the spectrum (the same frequency can be observed in the figure 5 a)). As well known (Castaing et al. 1989) this peak corresponds to the large scale circulation frequency (f LS ≈ V f f /4H), which can be estimated from the free fall velocity V f f = √ gβ∆T H which is ∼ 6 cm/s for lowest Ra H and ∼ 11 cm/s for the highest Ra H . In both single-phase and two-phase cases, a higher level of thermal power is seen for higher Ra H numbers. The same trend is seen for all the measurement positions at mid-height. If we now compare the single-phase and two-phase spectra, we can see that with the bubble injection, the thermal power of the fluctuations is increased by nearly three orders of magnitude. The bubbly flow also shows fluctuations at a range of time scales, with a gradual decay of thermal power from f 0.1 Hz − 3 Hz. The observation that substantial power of the temperature fluctuations resides at smaller time scales, as compared to the single-phase where the power mainly resides at the largest time scales, further confirms that the bubble-induced liquid fluctuations are the dominant contribution to the total heat transfer.
The thermal power spectra plots in figure 11 (a) show a Ra H dependence for both single-and two-phase cases. While upon normalising with the scale of temperature fluctuations T 2 rms in the single phase we do not observe complete overlap of the spectra possibly due to noise present at higher frequencies, in bubbly flow we observe a nearly perfect collapse of the three Rayleigh numbers (see figure 11(b)). This suggests a universal behavior for bubbly flow. Interestingly, this is similar to the velocity fluctuations spectra observed for bubbly flows (Lance & Bataille 1991;Riboux et al. 2010;Roghair et al. 2011;Prakash et al. 2016), where a normalisation with the scale of velocity fluctuations u 2 rms demonstrates universality. Furthermore, the same behavior is seen at all measurement positions and all Rayleigh numbers (note that in figure 11 (b), we have shown the measurements at the centre only). We also observe a clear slope of −1.4 at . Raw (a) and normalised (b) power spectra of the temperature fluctuations at the center of the setup for single-phase (red lines) and two-phase α = 0.9% (blue lines); solid line: RaH = 5.2 × 10 9 ; dashed line: RaH = 1.6 × 10 10 ; dash-dotted line: RaH = 2.2 × 10 10 .
the scales f 0.1 Hz − 3 Hz. It remains unclear why this exact slope is present, and how it can be attributed to bubble-induced turbulence. Events occurring at shorter time scales, such as at frequencies where the −3 slope is present in velocity spectra for bubbleinduced turbulence, typically starting at 1/t pseudo ∼ 35 Hz (Riboux et al. 2010), would be undetectable here due to the limiting response time of the thermistor used.

Summary of main results and discussion
An experimental study on heat transport in homogeneous bubbly flow has been conducted. The experiments are performed in a rectangular bubble column heated from one side and cooled from the other (see figure 1). Two parameters are varied: the gas volume fraction and the Rayleigh number. The gas volume fraction ranges from 0% to 5%, and the bubble diameters are around 2.5 mm. The Rayleigh number is in the range 4.0 × 10 9 − 1.2 × 10 11 . First, we focus on characterization of the global heat transfer for single-phase and twophase cases. We find that two completely different mechanisms govern the heat transport 10 10 10 11 10 12 10 13 10 14 Ra H 10 2 10 3 N u α = 0 % α = 0.5 % α = 5.0 % Figure 12. Extrapolation of the Nusselt number to high Rayleigh numbers for the two-phase (dashed lines -lowest and highest studied α), and the single phase case (dotted line). A crossover between the extrapolated single-and two-phase cases occurs at RaH ≈ 1.2 × 10 12 for α = 0.5%, and at RaH ≈ 4 × 10 13 for α = 5%. However, as we approach these RaH , we expect the Rayleigh-independent trends of Nusselt number to change, namely to increase with increasing RaH .
in these two cases. In the single-phase case, the vertical natural convection is driven solely by the imposed difference between the mean wall temperatures. In this configuration the temperature acts as an active scalar driving the flow. The Nusselt number increases with increasing Rayleigh number, and as expected effectively scales as: N u ∼ Ra 0.33 H (see figure  6 (a)). However, in the case of homogeneous bubbly flow the heat transfer comes from two different contributions: natural convection driven by the horizontal temperature gradient and the bubble induced diffusion, where the latter dominates. This is substantiated by our observations that the Nusselt number in bubbly flow is nearly independent of the Rayleigh number and depends solely on the gas volume fraction, evolving as: N u ∝ α 0.45 (see figure 8). We thus find nearly the same scaling as in the case of the mixing of a passive tracer in a homogeneous bubbly flow for a low gas volume fraction (Alméras et al. (2015)), which implies that the bubble-induced mixing is indeed limiting the efficiency of the heat transfer.
We further performed local temperature measurements at the mid-height of the setup for the gas volume fraction of α = 0% and α = 0.9%, and for Rayleigh numbers 5.2 × 10 9 , 1.6 × 10 10 and 2.2 × 10 10 . For single-phase flow, we observe that the mean temperature remains constant in the bulk at mid-height. However, in the two-phase flow case, this is completely obstructed by the mixing induced by bubbles. Injection of bubbles induces up to 200 times stronger temperature fluctuations (see figure 9 (b)). These fluctuations over a wide spectrum of frequencies (see figure 11) are thus the signature of the heat transport enhancement due to bubble injection. A clear slope of −1.4 at the scales f 0.1 Hz−3 Hz was also observed. In order to understand why this slope is present in that range it would be beneficial to perform local velocity measurements in the flow with heating, which is objective of our future studies. Further examination with fully resolved numerical simulations will also help us understand the effect of bubbles. These simulations are planned for future work, as well.
To conclude, we observe up to 20 times heat transfer enhancement due to bubble injection (see figure 6 (b)). This demonstrates that the diffusion induced by bubbles is a highly effective mechanism for heat transfer enhancement. Nevertheless, several questions remain unanswered. One question of great practical importance is: at what Rayleigh number will the contribution of natural convection to the total heat transfer become comparable or even greater than the contribution of bubble-induced turbulence? If we extrapolate the data to higher Rayleigh numbers, we can obtain the maximum expected value of Ra H at which the two contributions are comparable (see dashed lines in figure  12). This occurs at Ra H ≈ 1.2×10 12 for α = 0.5% and at Ra H ≈ 4×10 13 for α = 5%. As we approach these Ra H , we expect the Rayleigh-independent behavior of Nusselt number to change. Presumably, when Ra H is sufficiently large, the Nusselt number will increase with Ra H even for the two-phase cases. Based on our current knowledge, it is difficult to predict at what Ra H this trend will change. This calls for future investigations spanning a wider range of control parameters.