Modification of turbulence in Rayleigh–Bénard convection by phase change

Heavy or light particles introduced into a liquid trigger motion due to their buoyancy, with the potential to drive flow to a turbulent state. In the case of vapor bubbles present in a liquid near its boiling point, thermal coupling between the liquid and vapor can moderate this additional motion by reducing temperature gradients in the liquid. Whether the destabilizing mechanical feedback or stabilizing thermal feedback will dominate the system response depends on the number of bubbles present and the properties of the phase change. Here we study thermal convection with phase change in a cylindrical Rayleigh–Bénard cell to examine this competition. Using the Reynolds number of the flow as a signature of turbulence and the intensity of the flow, we show that in general the rising vapor bubbles destabilize the system and lead to higher velocities. The exception is a limited regime corresponding to phase change with a high latent heat of vaporization (corresponding to low Jakob number), where the vapor bubbles can eliminate the convective flow by smoothing temperature differences of the fluid.


Introduction
Multiphase fluid systems undergoing thermal convection are frequently encountered in industry and nature, from boilers and condensers to cloud and atmospheric dynamics. The prevalence of such systems has prompted interest in understanding the complex interaction between phase change and thermal convection and, particularly, how phase change affects the global properties of the flow. The standard and well-studied Rayleigh-Bénard (RB) cell [1]- [4] has been employed in recent numerical and experimental works to address questions about cloud formation in moist convection [5] and heat transport in the boiling process [6,7]. Experiments performed on ethane near its boiling point by Zhong et al [7] showed a significant enhancement of heat transport compared to single-phase transport, consistent with the numerical results from simulations of water near its boiling point performed by Oresta et al [6].
Here, we perform simulations of boiling in a cylindrical RB cell to gain further insight into how the phase change can modify the velocity and temperature fields, and turbulence level in thermal convection. As in single-phase RB convection, the dynamics are determined by the strength of the thermal forcing (the Rayleigh number) and the ratio of the kinematic viscosity to thermal diffusivity (the Prandtl number) [1]- [4]. The global response of the system is measured via the total heat transport through the cell (the Nusselt number, Nu) and the turbulence intensity (the Reynolds number, Re). For boiling, a critical additional parameter governing the heat transfer is the Jakob number, the ratio of the sensible heat to the latent heat of vaporization, Ja, which we vary to explore the different ways the phase change affects the response of the system.
When vapor bubbles form in a convecting liquid, it is not a priori clear how the velocity field and turbulence intensity will be modified. The dispersed bubbles have complex thermal and mechanical interactions with the liquid phase, see e.g. [8]- [17]. On the one hand, the density contrast between the liquid and vapor will induce motion due to buoyancy, but, on the other hand, the phase change from liquid to vapor removes energy from the liquid phase. It was proposed in [6], which focused on the physics of heat transfer in multiphase RB convection, that destabilization due to buoyancy dominates over stabilization due to thermal smoothing in most situations. In this paper, that idea is directly checked through calculation of the Reynolds number, with and without the thermal and mechanical feedback from the bubble on the flow.

The model and numerical method
Numerical simulations of water near its boiling point in an aspect ratio 1/2 (diameter/height) cylindrical cell are performed using the code developed by Verzicco and Camussi [18,19] and extended by Oresta et al [6] to include phase change. For this system, Pr = ν/κ = 1.75 and a temperature difference of T = 0.25 K is applied across a cell height of H = 17.9 mm, fixing Ra = gβ T H 3 /νκ = 2 × 10 5 (g is the gravitational acceleration, β the coefficient of thermal expansion of the liquid, ν the kinematic viscosity and κ the thermal diffusivity). The Jakob number, which determines how quickly bubbles will grow or shrink, is varied in the simulations and is Ja = ρc p (T h − T sat )/ρ V L (where c p is the liquid specific heat, T h the temperature of the hot bottom plate, ρ V the vapor density, L the latent heat and T sat = T h − T /2, the saturation temperature that is halfway between the hot and cold plate temperatures). For reference, water in these conditions would correspond to Ja = 0.37. The limit Ja = 0 implies an infinite latent heat so that bubbles cannot grow or shrink.
The code uses a finite difference scheme to directly solve the continuity, momentum and energy equations for the liquid phase in cylindrical coordinates on a non-uniform mesh, under the Boussinesq approximation [18,19]. Phase change is included by nucleating vapor bubbles randomly at the bottom plate so that the total number of bubbles in the cell, N , is fixed [6]. After testing the accuracy against results from finer grids, a resolution of 33 × 25 × 80 (in θ , r , z) is used. The bubbles move through the cell and grow or shrink depending on the local thermal interactions with the liquid. If a bubble condenses within the cell after encountering a cooler region, a new bubble is nucleated at the bottom plate at the next time step. If a bubble reaches the top plate without condensing, it is removed and a new bubble is nucleated, under the assumption that the timescale for condensation at the hot plate is very fast compared to the other system timescales 9 . An initial nucleation diameter of 25 µm is used, and the results for Ja = 0 do not depend strongly on this choice, due to the rapid growth in the hot thermal boundary layer at the bottom plate. In the limit Ja = 0, radial growth is restricted, so the initial diameter is important; however, the qualitative features described here are not affected by this choice, so the same nucleation diameter of 25 µm is used throughout this paper.
The bubbles are treated as point-like objects, and react back on the fluid as point sources of momentum and heat. Details of the numerical method and implementation can be found in [6] for the bubbles and in [18,19] for the Boussinesq equations; hence, only the essential equations are repeated here. The incompressibility condition (∇ · u = 0) is enforced on the liquid, and the momentum equation including the point forces from the bubbles is with u being the liquid velocity, p the pressure, T the variable liquid temperature and µ = ρν. The momentum forcing from the ith bubble at position ,i is the current radius [10]. A balance of added mass, drag, lift and buoyancy forces 4 determines the bubble motion [6,10,11]: where v i is the bubble velocity, and C A , C D and C L are the coefficients of added mass, drag and lift, respectively, the values of which are approximated as in [6]. The liquid energy equation including the contribution from the phase change is where Q i is the energy source or sink due to the ith bubble, which is proportional to the area of heat exchange, the temperature difference between the liquid and the bubble (assumed to be T sat ), and the heat transfer coefficient h b,i : The heat transfer coefficient depends on the single bubble Nusselt number, which is approximated as in [6]. In view of the small temperature differences, the bubble radial motion may be assumed to occur in conditions of quasi-equilibrium, so that the vapor pressure equals the ambient pressure and, therefore, the bubble temperature the saturation temperature. Thus, in this model, it is assumed that heat exchange goes completely into the phase change process, via the latent heat of vaporization, rather than increasing the temperature of the vapor. Hence, determines the rate of bubble growth.

Results
From the simulation results, we can directly calculate the Reynolds number of the flow from the volume-and time-averaged (denoted by brackets) root-mean-square (rms) velocity, u rms = ( u(r, θ, z) 2 V,t ) 1/2 . 10 The mean Reynolds number of the flow, calculated as Re = u rms H/ν, indicates the overall strength of the flow and the degree of turbulence in the system. In analyzing the flow structures, we will also look at the variation in the velocity field as a function of height using Re(z) based on u rms (z) = ( u(r, θ, z) 2 r,θ,t ) 1/2 . Simulations run until a statistically stationary state is reached, after which point the time-averaged quantities are computed over separate periods of time, with a variation in the quantities of less than 1% for the data in figure 1. Figure 1(a) shows Re as a function of N for Ja = 0.2. Introducing 1000 two-way coupled bubbles, corresponding to a total mean void fraction of only 7 × 10 −4 , leads to an increase of over an order of magnitude in the mean Re. Although these bubbles are injected with a diameter of 25 µm, they grow in the hot boundary layer and reach an average diameter of 115 µm. Increasing the number of bubbles causes further increase of the intensity of the flow, with more bubbles rising and driving the flow (figure 2(c)). However, the bubbles grow less as their number increases-for 10 000 bubbles, the mean diameter is reduced to 87 µm. This effect is likely due to a shorter residence time in the hot lower boundary layer (figure 4). In this way, the system is self-regulating, so that as more bubbles are introduced, the heat transfer and Re grow at a slowing pace. The initial bubble radius was not found to be an important factor in the global system properties, since smaller bubbles grow more near the hot boundary due to their smaller buoyancy and, once in the bulk, have the same size as those with larger initial radii.
Deactivating the mechanical (f i , equation (1)) and thermal (Q i , equation (3)) feedback from the bubbles separately reveals that the increase in Reynolds number is due to the mechanical forcing from the rising of the bubbles. Moderating the destabilizing mechanical effect is the thermal feedback, which stabilizes the flow and reduces Re. In fact, if the mechanical feedback is deactivated, Re quickly becomes O(10 −2 ) due to thermal smoothing of temperature gradients.
Smoothing of the temperature field occurs because the phase change tends to cool the hot regions of liquid, while it heats the cooler regions, weakening the thermal gradients that drive convection.
In contrast to the larger Ja = 0.2 vapor bubble case, if the bubbles are kept at a fixed size of 25 µm by setting Ja = 0 ( figure 1(b)), the mean Re decreases below the single-phase value. Fixing the bubble radius limits the mechanical feedback force f i , and also the contribution to the added mass force from the bubble's radial dynamics (proportional to dR b /dt) disappears [11].  In this case thermal smoothing dominates over the mechanical destabilization ( figure 2(b)). With few added bubbles, the thermal feedback stabilizes the flow, but as more bubbles are added flow is driven via buoyancy, and the Re begins to increase. It is possible to have a minimum, most stable configuration where the weakest flow occurs. This minimum near 5000 bubbles also corresponds to a flow configuration change from the typical horizontal roll to an axisymmetric toroidal roll (figure 2, also observed in [6]). Because the bubbles are introduced uniformly at the lower plate, when thermal smoothing wins the competition at Ja = 0, the system transitions to an axisymmetric state. Visualization of the stabilization via the thermal feedback is shown in figure 2 where snapshots of the dimensionless vertical velocity and temperature fields are shown in a vertical cross section through the axis of the cell. Both the reduction in velocity scales and the homogenization of the temperature field are clearly evident when 10 000 bubbles at Ja = 0 are introduced. At the other extreme, bubbles with Ja > 0 drive a strong convective flow due to their rising motion.
The vertical velocity variations as a function of height (plotted via the local Re z ) are shown in figure 3. Low-Ja-number bubbles have limited growth owing to the high latent heat of When Ja > 0, the bubbles grow rapidly near the hot bottom plate, gain speed and rise quickly through the top half of the cell. Some bubbles condense completely and disappear when encountering colder jets in the bulk region as well. This causes a noticeable asymmetry across the mid-plane, with the result that the vertical velocity variations are higher belowz = 0.5. The higher probability of finding a bubble near the bottom of the cell compensates for the fact that the rise speeds are lower there, and the non-dimensional bubble flux,q z =v z n, is still weighted towards the bottom of the cell ( figure 4(b)).

Conclusions
Phase change is demonstrated to have strong and varying effects on the velocity field and heat transport in thermal convection, depending on the ratio of latent heat to sensible heat (the Jakob number). By investigating separately the effect of the mechanical and thermal feedback of the bubbles on the fluid, we are able to directly see the competition between the two effects. In the low Ja limit, vapor bubbles extinguish motion and the cell becomes conductive, with a very low Reynolds number. For higher Ja, the mechanical forcing due to the rising motion of the bubbles drives additional velocity fluctuations in the liquid and leads to higher heat transport and increased Reynolds number of the flow. The asymmetric spatial variation of the fluid velocity field can be understood by looking at the bubble concentration and velocities-where the bubble vertical flux is highest, the vertical velocity variations in the flow are enhanced.
A straightforward nucleation method was used in this work, with bubbles nucleated at random positions on the hot plate. However, one can easily imagine that preferential nucleation with a non-uniform distribution, as might occur in experiments, could drive an interesting dynamical coupling between the flow and bubbles and is a potential outlook for this work. Calculations of the local Nusselt number u z T (z) could also help us to answer, in analogy to the question of whether the plumes or the background carry the heat in single-phase RB convection [20], the question of whether the bubbles or the liquid is mostly responsible for carrying the heat when there is a phase change.