Turbulent Rayleigh-B\'enard convection with bubbles attached to the plate

We numerically investigate turbulent Rayleigh-B\'enard convection with gas bubbles attached to the hot plate, mimicking a core feature in electrolysis, catalysis, or boiling. The existence of bubbles on the plate reduces the global heat transfer due to the much lower thermal conductivity of gases as compared to liquids and changes the structure of the boundary layers. The numerical simulations are performed in 3D at Prandtl number Pr=4.38 (water) and Rayleigh number $10^7\le Ra \le10^8$. For simplicity, we assume the bubbles to be equally-sized and having pinned contact lines. We vary the total gas-covered area fraction $0.18 \le S_0 \le 0.62$, the relative bubble height $0.02\le h/H \le0.05$ (where $H$ is the height of the Rayleigh-B\'enard cell), the bubble number $40 \le n \le 144$, and their spatial distribution. In all cases, asymmetric temperature profiles are observed, which we quantitatively explain based on the heat flux conservation at each horizontal section. We further propose the idea of using an equivalent single-phase setup to mimic the system with attached bubbles. Based on this equivalence, we can calculate the heat transfer. Without introducing any free parameter, the predictions for the Nusselt number, the upper and lower thermal boundary layer thicknesses, and the mean centre temperature well agree with the numerical results. Finally, our predictions also work for the cases with much larger Pr (e.g. $400$), which indicates that our results can also be applied to predict the mass transfer in water electrolysis with bubbles attached to the electrode surface or in catalysis.


Introduction
In wall-bounded buoyancy-driven turbulence, the boundary conditions play a crucial role in the flow structure and the global transport property of the system. Bubbles attached to the wall affect these boundary conditions. They often occur in various industrial applications. One example is water electrolysis, where bubbles are generated at the electrodes and can significantly reduce the global mass transport of the system by reducing the active electrode area (Vogt & Balzer 2005;Wang et al. 2014;Yang et al. 2018;Sepahi et al. 2022), leading to the decrease of the electrolyser efficiency. Another example is catalysis, where bubbles are generated by chemical reactions and can block the catalytic surface, thus also reducing the mass transport (Somorjai & Li 2010;Oehmichen et al. 2010;Xu et al. 2018). One example from daily life is heating water. When a pot of water is heated from below, many tiny gas bubbles nucleate at the bottom wall, reducing the heat transfer efficiency of the system due to the much lower thermal conductivity of gas as compared to the liquid. In all of these examples, bubbles attached to the wall influence the boundary layer (BL) and thus affect the global heat or mass transfer of the system. Therefore, it is highly desirable to quantitatively understand how much the global transport properties are changed. This is the motivation for our study.
As model system we choose Rayleigh-Bénard (RB) convection which is the paradigm of thermally driven turbulence (see the reviews of Ahlers et al. (2009) ;Lohse & Xia (2010); Chillà & Schumacher (2012); Shishkina (2021)), where a fluid between two parallel plates is heated from below and cooled from above. In previous studies, the effect of various plate properties on the flow structure and the heat transport were examined, e.g., staggered conducting and insulating strips on the plate (Wang et al. 2017;Bakhuis et al. 2018), temporally-modulated temperature (Jin & Xia 2008;Yang et al. 2020), plate with roughness (Zhu et al. 2017;Jiang et al. 2018;Zhu et al. 2019), and plates with different wettabilities (Liu et al. 2022). Here we pick RB convection as model system to study the effects of attaching bubbles on the global heat or mass transport in buoyancy-driven turbulence.
We employ direct numerical simulations with an advanced finite difference method combined with the phase field method (Liu et al. 2021). The bubbles are put at the lower (hot) plate with pinned contact lines so that they cannot move or detach from the plate. To focus on the effects of the Rayleigh number (dimensionless strength of the thermal driving) and the bubble geometry, we disregard mass exchange between the bubble and the liquid, and keep the bubbles at constant volume. The calculations are performed for various Rayleigh numbers and geometries, expressed through the relative area covered by the bubbles, bubble height, bubble number, and type of bubble distribution. To define the thermal BL thickness in the multiphase system, we extend the traditional temperature profile slope method, based on heat flux conservation. With this, we propose an equivalent single-phase system, for which we apply the Grossmann-Lohse (GL) theory (Grossmann & Lohse 2000, 2001Stevens et al. 2013) to predict the heat transfer and the temperature profile for RB convection with bubbles attached to the hot plate. Without introducing any new free parameter, these predictions well agree with our numerical results. They work for both moderate and large Pr , which indicates that they can be applied to predict both the mass transport in water electrolysis with bubbles attached to the electrode surface and to catalytic surfaces on which bubbles have formed.
The organization of this paper is as follows: The numerical method and setup are introduced in Section 2. The flow features and heat transfer are shown in Section 3. We define the thermal BL thicknesses of this new two-phase system in Section 4, and propose an equivalent single-phase system in Section 5 to calculate the heat transfer. The paper ends with conclusions and an outlook.

Numerical method and setup
The three-dimensional simulations are performed in a cubic domain of dimensions H 3 . The numerical method (Liu et al. 2021) combines the phase-field method (Jacqmin 1999;Ding et al. 2007;Liu & Ding 2015) and an advanced finite difference direct numerical simulation solver for the Navier-Stokes equations (Verzicco & Orlandi 1996;van der Poel et al. 2015), called AFiD. Numerical details, validation cases, and convergence tests were already presented in our previous study (Liu et al. 2021).
The phase field method is widely used in simulations of multiphase turbulent flows (Soligo et al. 2021), where the liquid-gas interface is represented by contours of the volume fraction C of the liquid. The corresponding volume fraction of gas is 1 − C.
The evolution of C is governed by the Cahn-Hilliard equation, where u is the flow velocity, and ψ = C 3 −1.5C 2 +0.5C −Cn 2 ∇ 2 C the chemical potential. We set the Péclet number Pe = 0.9/Cn and the Cahn number Cn = 0.75∆x/H with ∆x being the mesh size and H being the height of the RB cell. The parameters Pe and Cn are taken according to the sharp-interface approach proposed in Ding et al. (2007); Yue et al. (2010); Liu & Ding (2015). The flow is governed by the Navier-Stokes equation, the heat transfer equation, and the incompressibility condition, where θ is the dimensionless temperature, F st = 6 √ 2ψ∇C/(CnWe) the nondimensionalized surface force, where W e = ρ l U 2 H/σ is the Weber number, with the surface tension σ. The vector G = {[C + Λ β Λ ρ (1 − C)] θ −ρ/Fr} z represents the dimensionless gravity. All dimensionless material properties (indicated by a tilde,q) are defined in a uniform way,q = C + Λ q (1 − C), where Λ q = q g /q l is the ratio of the material properties of gas and liquid, marked by the subscripts g and l, respectively. The global dimensionless parameters controlling the flow are listed in Table 1. The most important response parameter of the system is the heat transfer, which is quantified by the Nusselt number Nu = Q/(k l ∆/H), with Q being the dimensional heat flux.
The values of the control parameters are chosen mainly based on the properties of air and water (see Table 1), though for better numerical efficiency we take the density ratio Λ ρ = 0.01, about 10 times larger than in reality. Note that the exact value of this parameter hardly affects our results. Since the bubbles are pinned and the buoyancy and surface tension forces applied on bubbles are always balanced, we take F r = 1 and W e = 100, also for numerical convenience. The geometrical parameters of the bubbles are the relative covered area 0.18 S 0 0.62, the non-dimensionalized height 0.02 h/H 0.05, their number 40 n 144, and the spatial distribution (uniform, random The temperature field is color-coded and the interface between liquid and gas is marked in green. As seen from the mean temperature profile in (c), due to the bubbles the mean centre temperatureθ = 0.461 is lower than the average temperature 1/2 of the top and bottom plates (dashed black line). and half-covered). Note that from S 0 ,h = h/H and n, we can also calculate the bubble volume (= S 0h /2+nπh 3 /6), the bubble contact radius (= S 0 /(nπ)), the bubble contact Since the local Weber number of the bubble is relatively low due to the small bubble size, the bubbles maintain their spherical-cap shape although the bubbles could deform and the flows inside and outside the bubbles are both solved. The bubble shape is closer to spherical with largerh and smaller S 0 /n.
To ensure that the contact lines are pinned on the plate, we set the value of C on the hot bottom plate equal to the initial value as the boundary condition of the phase field. The other boundary conditions are no-slip velocities on the top and bottom plates, fixed temperature θ cold = 0 (top) and θ hot = 1 (bottom), and periodic conditions in the horizontal directions. Stretched grids with 288 3 gridpoints are used for the velocity and temperature fields, and uniform grids with 576 3 gridpoints for the phase field (Liu et al. 2021), which corresponds to at least 12 gridpoints used for the bubble height. The mesh is sufficiently fine and is comparable to corresponding single-phase studies (Stevens et al. 2010;van der Poel et al. 2013).

Flow features and heat transfer
A typical thermal structure in RB convection with bubbles attached to the hot plate is shown in figure 1. Although there is continuous emission of thermal plumes, the bubbles almost maintain the shape of spherical cap due to the pinned contact lines and sufficiently strong surface tension. This is indeed the relevant situation during most of the time for the bubbles in electrolysis or catalysis. In figure 1(a) and (b), we observe the plumes rising up from the gaps between the bubbles. Inside the bubbles (see the inset in figure 1c), pure thermal conduction takes place since the local Rayleigh number inside the bubble Ra g ≈ β g gh 3 (∆/2)/(ν g κ g ) is small enough to remain under the onset of convection Ra c ≈ 1708 (namely Ra g < 100) for all cases. Furthermore, considering that the thermal conductivity of gas is much lower than that of liquid (i.e. for their ratio Λ k = 0.042 ≪ Rayleigh number Ra = β l gH 3 ∆/(ν l κ l ) 10 6 − 10 14 10 7 − 10 8 Prandtl number 1), the heat transfer through the gas phase is negligible. Therefore the overall heatconducting ability of the fluid near the hot bottom plate is lower than that near the cold top plate. Consequently, since the total heat flux is the same across each horizontal plane, the temperature drop across the hot bottom BL becomes larger than that across the cold top BL, leading to the mean centre temperature smaller than 0.5 (e.g. 0.461 in figure 1c), closer to the temperature of the cold top plate to that of the hot bottom plate.
In figure 2, we plot the heat transfer Nu (normalized by Nu s of single-phase system) as function of the geometrical parameters. As shown in figure 2(a), unsurprisingly, Nu/Nu s decreases with increasing the relative bubble-covered area S 0 due to the decreasing conducting area. We also observe that Nu/Nu s decreases with increasing bubble height h (normalized by the thermal BL thickness λ) in figure 2(b), since for larger bubbles there is less liquid (the major conducting fluid) in the thermal BL. With increasing h/λ, the simulation results gradually deviate from the predictions, since the larger h/λ, the more parts of the bubbles enter the bulk, which gradually deviates from our assumption that bubbles only affect the thermal BL structure. In contrast to the significant effects of S 0 andh on the heat transfer, Nu/Nu s is insensitive to the bubble number n (see figure 2c), since n only contributes to the second term in the bubble volume (= S 0h /2 + nπh 3 /2), which is of order O(h 3 ) and thus negligible compared to the first term (O(h)). We further show the effects of the spatial bubble distribution, including uniform, random, and the half bubble-covered distributions in figure 2(c). With all three types of distribution, the values of Nu/Nu s are almost the same. This indicates that the overall heat flux is insensitive to the spatial bubble distribution, at least in our model system. This differs from the previous studies (Liu & Zhang 2008;Jiang et al. 2018), where the large-scale flow is changed by the solid elements on the plate. The rigid solid surfaces imply no-slip velocity boundary conditions, which affect the flow structure much more than the boundary conditions of continuous velocity and shear stress on the deformable bubble interfaces.
Analogous trends are also found in the relationship between the mean centre temperatureθ and the geometrical parameters, as shown in figure 3. In all cases, the temperature The red lines denote the temperature profiles, the black lines the mean centre temperatureθ, the blue lines the solution for pure thermal conduction in the BL, i.e. eq. (4.1), and the green dashed line the tangent of (4.1) at the intersection. θ * is the intercept of the tangent. We define the dimensional BL thickness near the hot and cold plates as λ hot and λ cold , respectively, i.e., the distance from the black dashed lines to the corresponding plate.
profile is asymmetric due to the bubbles, namelyθ is lower than the average temperature 1/2 between top and bottom plates, as explained above. The temperature profile can be quantitatively described byθ and the top and bottom thermal BL thicknesses, the values of which are calculated in Section 5. Such an asymmetric temperature profile has already been studied in the context of single-phase RB convection under non-Oberbeck-Boussinesq (NOB) conditions (Ahlers et al. 2006), where in water the viscosity and thermal diffusivity are temperature dependent and smaller at the hot bottom plate than at the cold top plate. Also this leads to an asymmetric temperature profile. Ahlers et al. (2006) employed an extended Prandtl-Blasius BL theory for the NOB conditions in the two BLs, and coupled them by imposing heat flux conservation. With this they could calculate the centre temperature and the thermal BL thicknesses, which well agree with the experimental results. Note that that analytical calculation must be adapted to be applicable to the situation here, since here the BL flow is not parallel to the plates due to bubbles. In addition, whereas in the NOB case one has no-slip and no-penetration boundary conditions throughout, here on the spherical-cap shaped bubble interface we have the emerging condition of the velocity and shear stresses continuity. However, as in Ahlers et al. (2006), we use the concept of heat flux conservation to define the thermal BL thickness (see section 4), and then propose the idea of an equivalent single-phase system to mimic the system with attached bubbles (see section 5).

Thermal BL thicknesses
A sketch of the temperature profiles near the plates is displayed in figure 4, where we also show how the thermal BL thicknesses are defined. Near the cold top plate (see figure  4a), the thermal BL thickness λ cold is defined through the usual convenient definition via the slope of the temperature profile at the plate from assuming a pure conductive thermal BL and a well-mixed bulk (Ahlers et al. 2006). As λ cold we take that distance from the plate, where the tangent to the temperature profile at the plate reaches the mean centre temperatureθ. However, this definition cannot be directly applied near the hot bottom plate, due to the bubbles. As the bubbles reduce the area occupied by the liquid, the conducting area varies with z. We therefore correct the solution for pure thermal Ra * = 2θRa θ sim Figure 6. (a) Comparison between the mean centre temperatureθsim and the effective temperature θ * sim , both taken from simulations. The symbols are the same as in figure 2. The shadow is the zone within ±5%. (b) Sketch of the bubble-system under consideration and the equivalent single-phase RB system.θsim and θ * sim are plotted this way to support the existence of the equivalent single-phase system. conduction in the hot bottom BL based on heat flux Q conservation, Here, the term on the right is the heat flux on the cold top plate. The gas-covered area S(z, S 0 ,h, n) = max[(S 0 + nπzh)(1 − z/h), 0] is the insulating area, which is defined with the assumption of spherical-cap shaped bubbles. The thermal BL thickness near the hot bottom plate λ hot equals the distance from the plate to the intersection between the lines of (4.1) and of θ =θ, as shown in figure 4(b). With the definitions above, both values of λ hot and λ cold are close to each other in all cases (within 10%), as shown in figure 5. Since λ hot and λ cold are almost the same and θ < (θ hot + θ cold )/2, the temperature variations across the hot bottom BL and the cold top BL are different, which is reflected in that the temperature profile is asymmetric.

Predictions for
Nu and the centre temperature using equivalent single-phase RB system To calculate the heat transfer Nu and the mean centre temperatureθ, we propose the idea of using an equivalent single-phase setup to mimic the system with attached bubbles. To find the equivalent flow, we first obtain the effective temperature at the hot bottom plate. This is done by plotting the tangent of (4.1) at the position of λ hot , and the obtaining θ-intercept of the tangent as effective temperature θ * (see figure 4b).
Next, we compare the effective temperature θ * and the mean centre temperatureθ in figure 6(a), which shows that θ * approximately equals to 2θ (within ±5%). Thus we assume a nearly equivalent single-phase RB counterpart with θ hot = 2θ and θ cold = 0, such thatθ is the mean value of temperatures on the two plates, as shown in figure 6(b). Again, we note that λ hot is close to λ cold , based on our definition of the thermal BL thickness in Section 3. Thus, we have the following relations: Here λ s is the thermal BL thickness in the single-phase system, which can be calculated as where Nu e is the heat transfer estimated from the GL theory (Grossmann & Lohse 2000, 2001 for the equivalent single-phase system with Ra * = 2θRa, θ hot = 2θ and θ cold = 0. In the two-phase system and the equivalent single-phase system, the dimensional heat transfer Q (= Nu k(θ hot − θ cold )/H) should be the same, which yields Nu(Ra, Pr , S 0 ,h, n) = 2θNu e (Ra * , Pr ), where Nu(Ra, Pr , S 0 ,h, n) = (θH)/λ cold is the heat transfer for the two-phase system. Combining the relations (5.1)−(5.3) for the equivalent single-phase RB system with the relationship (4.1) betweenθ and λ cold in section 4, we can now calculate the heat transfer in the system with attached bubbles.
We emphasize that with this approach, using the equations above without introducing any free parameter, for given Ra, Pr and bubble geometries, we can now calculate Nu, θ, λ hot and λ cold . The good agreements between the simulations and the predictions for Nu,θ, λ hot , and λ cold are shown in figures 2, 3, and 5, respectively, where all deviations between simulations and predictions are within ±5%.
We further check whether our approach is also applicable for large Pr , which, as explained above, has relevance to transport phenomena in water electrolysis and catalysis. Water electrolysis and catalysis can both lead to natural convection driven by buoyancy, which originates from the density difference of the solute with different concentrations of the electrolysis or catalysis product. Here, the mass transfer is also characterized by Nu, i.e. the mass transfer normalized by that with pure diffusion (in this context normally called Sherwood number Sh). The control parameters are the Grashof number (dimensionless strength of the solute driving) Gr and the Schmidt number (the ratio of viscous diffusion and mass diffusion rates) Sc, corresponding to Ra/Pr and Pr in RB convection, respectively. The value of Sc in water electrolysis is always large, e.g., Sc = 400 (Sepahi et al. 2022).
We performed two simulations (for two different bubble heights) at Pr = 400 and Ra = 10 7 (tabulated in Table 2) with a sufficiently fine mesh as explained in section 2. Note the much higher computational costs at this large Pr , due to the required long time (∼ Pr 1/2 ) for the system to enter the statistical steady state. The bubble geometries are characterized byh = 0.03 and 0.05, S 0 = 0.32, and n = 40. The good agreements between our parameter free predictions and the results from the simulations are shown in Table 2. This supports that our predictions can be directly applied to water electrolysis and catalysis.

Conclusions and outlook
Turbulent RB convection with gas bubbles attached to the hot plate is numerically investigated for 10 7 Ra 10 8 and Pr = 4.38 and 400. The bubble geometrical parameters are the relative bubble-covered area S 0 , the relative bubble heighth, the bubble number n, and the spatial bubble distribution. Due to the much lower thermal conductivity of gas as compared to liquid, the temperature profile is asymmetric and the heat transfer efficiency of the system is reduced. More specifically, Nu significantly decreases with increasing S 0 andh, but is almost unaffected by n and the types of bubble distribution.
To predict the heat transfer and the mean centre temperature of the system, we have proposed the idea of using an equivalent single-phase system to mimic the system with attached bubbles. By applying the GL theory for the equivalent system and imposing heat flux conservation in the two thermal BLs, we can predict the heat transfer, the top and bottom thermal BL thicknesses, and the mean centre temperature, without introducing any free parameter. The predictions well agree with the results from the simulations. Briefly, in Section 4 we got one relationship betweenθ and λ in eq. (4.1), and in Section 5 we got another relationship betweenθ and λ in eqs.(5.1) and (5.3). Then, for only given Ra, Pr and bubble geometries, we can well predict the heat transfer and temperature profile in the system with bubbles attached to the bottom plate.
The results of this study can be used not only for the heat transfer in RB convection with bubbles attached to the plates, but also for the mass transfer in electrolysis or catalysis. Our predictions can help to obtain estimates for relevant applications and e.g. optimize the heat or mass transfer and flow features in systems in which bubbles are forming on the plate(s). It would also be interesting to extend our basic idea to other wall-bounded turbulent systems with various plate properties, such as plates with inhomogeneous properties (e.g. wettability or conducting ability).