Wind Reversal in Bubbly Natural Convection

: The multi-phase Rayleigh–Bènard convection has been weakly investigated, even though it plays a leading role in the theoretical and applied physics of the heat transfer enhancement. For the case of moderate turbulent convection, a rather unexpected result is an unusual kind of wind reversal, in the sense that the ﬂuid is found to be strongly inﬂuenced by the bubbles, whereas the bubbles themselves appear to be little affected by the ﬂuid, despite the relative smallness of the Stokes numbers. The wind reversal induced by the bubbles dispersed in the ﬂuid is a new and remarkable phenomenon in multi-phase ﬂows that provides further perspectives in understanding the complex physics leading the enhancement of thermal convection. For this reason, the fundamental research proposed in this paper aimed to identify a space of control parameters and the physical mechanisms responsible for the wind reversal induced by dispersed bubbles in a conﬁned convective ﬂow. The strength of the following description lies in an innovative numerical approach, based on the multi-scale physics induced by the coupling of the local thermal and mechanical mechanisms arising between each bubble and the surrounding ﬂuid. The continuous phase has been solved numerically using the direct numerical simulation (DNS) technique and each bubble has been tracked by means of a particle Lagrangian model.


Introduction
The standard Rayleigh-Bénard convection occurs in a Newtonian single-phase fluid, bounded by lower (hot) and upper (cold) horizontal plates both at constant temperature, and an adiabatic side-wall. The thermal-induced buoyancy force pushes the stream toward the hotter region where the fluid is lighter. In the colder region, it happens in exactly the opposite way. The thermal convection undergoes significant changes with phase transitions [1,2] and when polymers [3], bubbles [4,5], and particles [6] disperse in the continuous phase. The most relevant effect is the heat transfer enhancement of the multi-phase convection respect to the case of single phase convection.
The multi-phase Rayleigh-Bénard convection plays a leading role in a wide range of physical phenomena such as the formation of atmospheric precipitation [7], magma chambers [8,9] boiling of liquid [10], and counterflow cooling towers [11].
It is well known that suspended particles or bubbles have a very strong effect on turbulent flows. A considerable amount of experimental evidence has been accumulating for more than 40 years, showing that small particles or bubbles dampen the turbulence, while larger ones increase it, due to vortex shedding behind them [12]. The particles selectively interact with the coherent structures and, in turbulent flows, they preferentially distribute avoiding strong vortical regions whereas they clustering T , is a negligible fraction of that at wall and essentially all the T are produced within the thermal boundary layers.
The strength of the numerical approach is the multi-scale physics induced by the coupling of local thermal and mechanical interactions between each particle and the surrounding fluid. For this purpose, the standard point-particle model has been extended in order to interact with the thermal effects associated with the thermal balance at the particle surface. A significant difference respect to the single-phase numerical method lies in the new source terms introduced in the momentum and in the thermal energy equations. The time integration of the Navier-Stokes equations with the Boussinesq approximation and the Lagrangian particle equations takes into account the Kolmogorov time scales, the particle relaxation time for the momentum and the thermal balance. A new feature is a further constraint on the time step since the particle moves through the viscous and the thermal boundary layers with a sequence of thermal and mechanical equilibrium states with the surrounding fluid. The numerical results are in good agreement with experiments in the case of bubbly convection reported in [26]. The solver handles the multi-dispersed particle distribution and it is parallelized by OpenMp for both the Eulerian continuum phase and the Lagrangian particle tracking.

Methods
In a cylindrical domain, 7.5 × 10 3 bubbles were inserted so that each of them gets a random value for its initial position. The cylindrical domain is 17.9 mm height (H) with diameter of D = 8.95 mm and aspect ratio Γ = D/H = 0.5. The cold temperature at the lower plate is T c and the hot temperature at the upper plate is T h , which is ∆T = T h − T c = 0.25 K, whereas the dimensionless temperature is ψ = (T − T c )/(T h − T c ). The characteristic free fall velocity is U f = gβ(T h − T c )H and it is used to adimensionalize the velocities. The length scale is the height of the domain H and the time scale is τ = H/U f . The number of bubbles is constant and it is equal to N b = 10 5 . The bubbles change their volume because of the hot temperature and the maximum void fraction is equal to 5 × 10 −3 .
The main control parameters for the bubbles size, whose initial diameter is equal to 25 µm, are the Stokes number St f = 2.94 × 10 −5 , the Froud number Fr = √ 1.88 × 10 −2 , and the Jakob number Ja = 1.2 × 10 −1 . They disperse in a stationary convective flow with Rayleigh number Ra = 2 × 10 5 and Prandtl number Pr = 1.65 (water at a temperature of 373 K and atmospheric pressure) in order to investigate the change of the flow without background turbulence but subjected to a thermal gradient.

Continuum Phase
We solve the continuity, momentum, and energy equations in a cylindrical domain with a No-slip condition applied to the all solid surfaces. The temperatures of the top and bottom surfaces are kept constant at T c and T h , respectively, while the lateral boundary is assumed to be adiabatic.
The Navier-Stokes equations are solved in a cylindrical coordinate system using a second-order, finite difference, fractional-step method on a staggered grid.
in which u is the liquid velocity field. The momentum equations is where D/Dt is the convective derivative, p and T are the pressure and temperature, and µ = νρ the dynamic viscosity. The continuity equation reads as in Equation (1) in the case of small bubble volume fraction. The Navier-Stokes equations are solved using a finite difference numerical code in cylindrical coordinates based on direct numerical simulation (DNS). Direct numerical simulation is a methodology in computational fluid dynamics in which the whole range of spatial and temporal scales of the coherent structures are resolved without any turbulence model. The flow field is calculated by means of the fractional step method and an approximate-factorization technique, in which the non-linear terms are treated explicitly, whereas the viscous terms are computed implicitly. The fractional step is a method of approximation of the evolution equations based on the decomposition of the operators. The role played by the pressure is a projection operator that projects an arbitrary velocity field into a divergence-free velocity field. The general procedure for projection methods is a predictor-corrector approach. In the first step, a preliminary velocity field is computed using the momentum equation. This velocity does not satisfy the continuity equation because in this step the pressure value is not updated. In a second step, a Poisson-type equation for the pressure is solved, which is derived using the continuity equation. In the last step, the preliminary velocity field is projected in a divergence-free velocity field using the computed pressure [27]. We used the second-order explicit Adams-Bashforth scheme for the convective terms and the second-order-implicit Crank-Nicholson scheme for the viscous terms. Implicit treatment of the viscous terms eliminates the numerical viscous stability constraint. This constraint is strong for the flow with low-Reynolds number and it is relevant near boundaries where stretched meshes are used. The implicit treatment of the viscous terms requires the inversion of large sparse matrices. In order to face this numerical problem, the viscous terms are split into three tridiagonal matrices by a factorized procedure. The time advancement of the solution is obtained by a 3rd order Runge-Kutta. The simulations were carried out using a grid with 33 × 25 × 80 nodes in the azimuthal, radial and axial directions, respectively. We checked the grid resolution with a finer grid of 33 × 40 × 120 nodes finding differences in the value of the Nusselt number within 1.5%, [28]. Moreover, the viscous and the thermal layers were resolved since the points clustered toward each boundaries are 15 in the case of single phase flow and 6 in the case of multi-phase flow.

Dispersed Phase
The bubbles were considered to be point-like. The single bubble is located at the coordinate x i and the force f i that it induce on the continuum phase reads as follows: in which R bi is the radius of the i-th bubble. The continuum phase energy equation is: where k = κρc p is the liquid thermal conductivity and Q i is the energy source or sink due to phase change of the i-th bubble. We take into account the thermal balance between each bubble and the surrounding fluid by means of a heat transfer coefficient h bi where is the liquid temperature evaluated at the position of each bubble center. In the previous equation we assume that the phase change is slow and the temperature at the bubble surface is kept at the saturated conditions.
The equation of motion for each bubble balances added mass, drag, lift and buoyancy, in which C A , C L , and C D are the added mass, lift, and drag coefficients, respectively. The bubble radius R b is calculated by balancing the latent heat associated to evaporation or condensation with the heat exchanged with the liquid We take C A = 1/2, the standard potential-flow value for a sphere and we model the drag coefficient as suggested by [29]: The equation of the dispersed phase and continuum phase are solved numerically according to the following procedure: Upon integration of the fluid momentum and energy equations over a computational cell, the effect of the particles is localized at their position and this effect must be replaced by an equivalent one located at the computational nodes. This task is achieved by a second-order-accurate interpolation. The interpolation preserves the resultant and the couple of the particle forces, as well as the total amount of heat that each particle exchanges with the fluid. The numerical integration of the dispersed phase follows three steps: (i) the equation of the thermal balance at the vapor surface is solved by means the 2nd order Adam-Bashfort scheme and the radius of the bubble is updated; (ii) the momentum equation is solved using the 2nd order Adam-Bashfort scheme too, whereas the implicit 2nd order Crank-Nicholson method has been used for the drag force (iii) the bubble position is updated by Adam-Bashfort 2nd order. The implicit 2nd order scheme for the bubble momentum equation makes the solution stable with the smallest bubbles since they have a fast time evolution.
When a bubble reaches the top of the cell, it is removed from the calculation and a new bubble is re-injected at a random position on the bottom plate with a velocity equal to that of the surrounding fluid.
Each bubble is injected with its initial diameter as it is not estimated using any nucleation model. The mesh independence study and both the experimental and numerical validations were reported in Lakkaraju et al. [26]. More details were reported by the authors-in the numerical code by Oresta et al. [28].

Convection Enhancement Induced by Vapor Bubbles Nucleation
The control parameter quantifying the convection is the time evolution of the Nusselt number Nu (heat flux transferred through the plates normalized by its purely conductive value) as shown in  The Jakob number is a dimensionless parameter that takes into account the bubble ability to increase its own volume due to the heat transfer with the surrounding fluid. To increase the Jakob number, the bubbles enlarge more and more, traveling in the vertical direction towards the upper plate. Under this condition, the presence of bubbles increases the time averaged Nusselt number from 4.718 (no bubble flow) to 13.13 (two-way coupling flow). It should be pointed out that the bubble becoming bigger compared to its initial size changes its dynamics, which can be described by means of the point-like approximation. Moreover, the bubbles are responsible of the fluctuations in time of the Nusselt number, as their irregular rising motion generates small vortical structures characterized by highly non-stationary regime with high frequencies of the control parameter signal.
The dispersion of bubbles causes a strongly non-linear change in the heat transfer. In Figure 1 two flow cases, with and without bubbles, are shown for a computational domain with aspect ratio Γ = 1/2. The average value of the heat transfer trend of the flow with bubbles is much enhanced, with respect to non bubble flow and it is worth noting how the fluctuations determined the transition to a turbulent-like behavior. In Figure 1b, the probability density function of the Nusselt number fluctuations is shown As we see, the probability density function of Nu nicely agrees with the Gaussian distribution. This statement underline the role played by the rare events of the distribution tail superimposed to the mainstream circulation. The rare events due to the nucleation, growth, and rising motion of the bubbles, are independent because ruled by the Gaussian distribution. The bubbles are not clustered by the flow if we consider the average time at which any random configuration is allowed. Instead, if we look at the instantaneous displacement of the bubbles (Figure 2) they are preferentially located into the ascending stream. How can the bubbles behave differently? The answer is in the wind reversal, which instantaneously pushes the particles to the ascending stream, but also sweeps over the azimuthal direction, and in doing so, avoiding any preferential marker in the average thermodynamical parameters. In fact, they follow the Gaussian distribution such as the Nusselt number.
The following investigation was based on the evidence that the Nusselt number shows higher frequencies mainly because of small structures induced by the bubbles rather than the large-scale circulation.

Wind Reversal: Detection and Analysis
The bulk dynamics is characterized by the angular orientation θ c of the large-scale circulation (LSC) plane and by the amplitude of the temperature fluctuation δ in the azimuthal direction. Typically, δ is a measure of the magnitude of the LSC.
Taking advantage of the data accessibility provided by the direct numerical simulation, the flow dynamics were explored.
Several numerical probes are placed within the fluid volume, allowing point-wise velocity and temperature time series to be extracted. The eight probes (Figure 3) are located along the azimuthal direction every π/4 rad starting from θ 0 = 0 rad on the horizontal plane midway between the two plates (cold and hot) at a radial location r/H = 0.125 (where, r the radial coordinate and H the height of the cylindrical domain). According to Brown et al. [30], we fitted the adimensional probe temperature ψ i i = 0, ...7, separately at each time, using the following function: where ψ m = 1/8 ∑ 7 i=0 ψ i is the average temperature of the eight probes. The δ and θ c are estimated by enforcing the minimum of the square error function In order to avoid the non-linearity of the error function due to the cosine of θ c (coming from Equation (9)), that equation is rewritten as follows.
(11) Figure 4 shows the ψ i and the corresponding fit ψ * i at three time instants with Figure 4a corresponding to the beginning of the phenomenon. The parameters ψ m , δ and θ c have been listed in Table 1.
The fitted temperature ψ * i shows a good agreement with the sampled temperature ψ i sampled; moreover, we observe that δ and θ c are strongly dependent on time. The reorientation of the LSC was observed in turbulent convection but its dynamics is still only partially understood. An interesting aspect of the LSC observed in literature is its spontaneous changing in direction, due either to an azimuthal rotation of the entire structure without a significant change of the flow speed, or the reverse flow direction without a significant rotation of the circulation plane. We are referring to the latter process as motion cessation since it determines a momentary vanishing of the flow speed, and to both combined processes as reorientation. The reorientation cause by both mechanisms corresponds to a continuous range ∆θ c .
The reorientation of LSC happens several times in the thermal flow forced by bubbles. We analyzed the time evolution of δ and θ c in a time interval close to 528.78 ( Figure 5). For this time instant, δ shows a minimum value corresponding to a quick variation of θ c . This means that the fluctuation of the temperature vanishing and simultaneously the LSC is reoriented. Another proof of this behavior is provided by the time evolution of the flow axial velocity recorded in two symmetric positions (probes 3 and 7): its sign changes at the instant for which δ shows a minimum in connection with the sharp variation of θ c . The presence of large-scale structures and their reorientation it is very clear if we observe the flow axial velocity plotted in the horizontal section of the cylindrical domain at z = 0.5 ( Figure 6). In the same figure, the orientation of the LCS is reported at time 521.8 and 530. The presence of the large-scale and its orientation suggest a strong relation between the effects of the bubbles and the bulk dynamics of the thermal flow at high Rayleigh numbers without bubbles.
The bubble distribution is highly inhomogeneous (Figure 2) showing clustering and preferential sampling of the up-flow region (Figure 2b) before the reversal involving a random re-mixing process of the bubbles. This phenomenon is confirmed by the strong variation of the Stokes number related to the growing and shrinking process of each bubble as it shown in Figure 7. In the same figure, the volume fraction has a small fluctuation with an averaged value of 0.2.

Description of the Flow Coherent Structures
A quantitative way of characterizing the mean flow is determining its energy content: by performing a Fourier transform of the velocity field, its spectrum of n azimuthal modes can be obtained as well as their energy contents; namely, if u(θ, r, z) is the velocity field, its FFT in the azimuthal direction return the matrix û(n, r, z), where n the azimuthal wavenumber; the azimuthal energy modes are then obtained by integrating û(n, r, z)û * (n, r, z) in r and z, for each wavenumber n, û * being the complex conjugate of û. This decomposition is particularly significant since the n = 0 mode corresponds to the axial symmetric toroidal structures, the n = 1 mode contains the energy of the large-scale structures spanning the cell, and the n ≥ 2 mode are the structures with n − f old azimuthal symmetry.
In Figure 8 time history of the ratio between the first and zero mode azimuthal energy: solid line indicates the no bubble flow whereas the dashed line refers to the two-way coupling case. The first mode has higher amount of energy respect to the zero mode only for the flow without bubbles. On the contrary, the dispersion of bubbles makes the amount of energy E 1 and E 0 comparable at every instant of time in the other case. This means that the large-scale is characterized by a big roll filling the whole domain that temporarily in time switches to structures characterized by a roll coexisting with an axisymmetric vortex. For instance, the flow exhibits these two flow structures at the time 530 when the E 1 /E 0 2.5 ( Figure 9). The vertical plane crossing the axis at azimuthal position θ c keep the vortex filling the whole domain: that corresponds to the first mode. Instead, the vertical plane Instead, the vertical plane crossing axis at azimuthal position θ c + π/4 keeps the axisymmetric vortex, and that corresponds to the zero mode.  In addition to the effect of the bubbles on the velocity field, the temperature field also deeply changes in the force of the moving bubbles, as shown in Figure 10. In the case of stationary thermal flow without bubble (Figure 10a), the thermal field shows a fluctuation of the temperature in the azimuthal direction that decreases with the two-way coupling because of its re-mixing. Moreover, the bubbles generate elongated structures attached to the upper plate commonly called "thermal plume". These plumes are like emissions of hot (cold) fluid from bottom (top) thermal boundary layers and are carried by their buoyancy. This observation is another aspect concerning the connection of the thermal flow at high Rayleigh numbers with the effects due to the bubbles dispersed in a stationary convective flow.

Conclusions
The Multi-phase Rayleigh-Bénard convection was investigated in a confined slender cylindrical domain. At moderate Rayleigh number, the single-phase flow shows a stationary state of the large-scale circulation characterized by a single roll filling the domain. When this coherent structure becomes unsteady, it causes the nucleation of vapor bubbles.
The time history of the Nusselt number shows a Gaussian distribution where the rare tail events connected with the bubble dynamics overlap with the average main stream flow.
This phenomenon apparently disagrees with the instantaneous clustering of the bubbles that preferentially place themselves in the ascending flow. Those effects coexist, since the LSC changes its azimuthal position in time, so sweeping the whole section area trough quasi-stationary azimuthal configurations.
The numerical approach could detect the wind reversal induced by the bubble dispersion showing that, by considering an averaged dynamics, the Nusselt number is randomly affected by the bubbles, even though it shows instantaneous preferential interaction with the flow structure populating the cylindrical domain.

Conflicts of Interest:
The authors declare no conflict of interest.