Time-Periodic Cooling of Rayleigh–Bénard Convection

: The problem of Rayleigh–Bénard’s natural convection subjected to a temporally periodic cooling condition is solved numerically by the Lattice Boltzmann method with multiple relaxation time (LBM-MRT). The study finds its interest in the field of thermal comfort where current knowledge has gaps in the fundamental phenomena requiring their exploration. The Boussinesq approximation is considered in the resolution of the physical problem studied for a Rayleigh number taken in the range 10 3 ≤ Ra ≤ 10 6 with a Prandtl number equal to 0.71 (air as working fluid). The physical phenomenon is also controlled by the amplitude of periodic cooling where, for small values of the latter, the results obtained follow a periodic evolution around an average corresponding to the formulation at a constant cold temperature. When the heating amplitude increases, the physical phenomenon is disturbed, the stream functions become mainly multicellular and an aperiodic evolution is obtained for the heat transfer illustrated by the average Nusselt number.


Introduction
The building sector is one of the largest energy consuming sectors. It represents a large proportion of total energy consumption, much higher than industry and transport in many countries [1][2][3][4][5]. The major challenge to guarantee good comfort is that energy consumption and the level of comfort are often in conflict in a room [6].
In a confined space, thermal comfort depends on the operation of controlled installations (i.e., heating, ventilation and air conditioning) [7]. Therefore, to assess the energy performance of the building, a thermal analysis is essential in order to predict the thermal responses and calculate building loads (heating/cooling). The key of good thermal comfort is the temperature control. To assess it, two models exist: physical or white-box models based on the energy and mass balance equations, and data-driven or black box models based on artificial neural network and developed after sufficient data are available.
In order to propose an adequate design of the indoor environment where people work and live, several papers have been published on the thermal comfort of buildings [8][9][10]. Rayleigh-Bénard (RB) convection is a classical problem of natural convection. It appears in several applications of engineering and building, such as thermal comfort by using floor heating or ceiling cooling. The other challenge is the energy storing, so either phase change materials or liquid storing are integrated to enhance the thermal building inertia.
Chandrasekhar [11] and Drazin and Reid [12] implemented a full report on the linearization theory. Moreover, many numerical studies of RB natural convection in rectangular enclosures have been carried out for Newtonian fluids [13][14][15]. Other researchers were interested in the Rayleigh-Bénard convection in viscoplastic fluids leading to numerical and experimental investigations [14][15][16][17]. The Rayleigh number (Ra) is the parameter that quantifies the intensity of the thermal driving in convection. For sufficiently large Ra, Rayleigh-Bénard convection flow becomes turbulent. Significant progress in our understanding of turbulent convection has been obtained by both experimental and numerical studies [18][19][20]. The Rayleigh-Bénard convection with all the issues illustrated above becomes more complex in the case of time-dependent boundary conditions [21]. The Rayleigh-Bénard problem, with linear temperature increase, was studied by Kaviany [22] and extended later on by Kaviany and Vogel [23], with the inclusion of solute concentration gradients. The solute gradient represents either the phase diagram solute redistribution near phase change interface or the stabilized long-term energy storing solar pound. According to the results of the Rayleigh-Bénard convection studies, a sufficient condition has been found to control the frequency of heat pulsation in order to initiate convection in a periodically heated and cooled cavity from top wall [24].
Aniss et al. [25] have studied the influence of the gravitational modulation on the stability threshold in the case of a Newtonian fluid confined in Hele-Shaw cell and subjected to vertical periodic motion. A time-dependent perturbation expressed in Fourier series has been applied to the wall temperature according to Bhadauria and Bhatia [26]. They found that it is possible to advance or delay the onset of convection. Umavathi [27] also investigated the effect of external modulation on the thermal convection in a porous medium saturated by a nanofluid. It was found in this last study that the low frequency symmetric thermal modulation is destabilizing while moderate and high frequency symmetric modulation is always stabilizing. Recently, Himraneet et al. [28] have studied numerically the Rayleigh-Bénard convection using the Lattice Boltzmann method. The authors considered periodic heating at the lower wall of the cavity. For high values of Rayleigh number, they obtained an unsteady regime in the form of temporal evolution with several frequencies. Abourida et al. [29] have examined Rayleigh-Bénard convection in a square enclosure with a top wall submitted to constant or sinusoidal cold temperature and sinusoidally heated bottom wall. It was observed that by varying the two imposed temperatures, basic differences were noted in comparison to the case of variable hot temperature and that of constant boundary temperature conditions. Raji et al. [30] have presented Rayleigh-Bénard convection inside square enclosure with a time-periodic cold temperature in the top wall. They proved that the variable cooling can lead to a significant improvement in heat transfer compared to constant cooling, particularly at certain low periods. The influence of thermoelectric effect on the Rayleigh-Bénard instability has been well investigated by researchers [31][32][33].
The literature review has shown that Rayleigh-Bénard's convection problem is still relevant and often considered with constant boundary conditions whereas reality suggests that heating and cooling conditions are modulated over time (i.e., scrolling of the days). In the present study, the problem of time-periodic cooling applied to Rayleigh-Bénard convection is investigated using Lattice Boltzmann method with multiple relaxation time (LBM-MRT). The study finds its interest in the field of thermal comfort where current knowledge has gaps in the fundamental phenomena (theoretical funds) requiring their exploration.

Mathematical Formulation
Lattice Boltzmann method was adopted for the resolution of fluid flow and transports phenomena. This method has non-linearity in its mathematical formulation and therefore approximates the temporal variation wall temperature (totally explicitly unsteady). Fluid flow is described as a movement of particles. In two-dimensional domain, the particles' moving is done by distribution functions with DnQm model, where "n" denote dimensions and "m" discrete velocities [34]. In our study, we used D2Q9 model for dynamic field and D2Q5 for thermal model. Finally, collisions and advection terms are computed by the Boltzmann equation of the distribution functions [35,36]: where f i is the distribution function with velocity e i at lattice node x at time t, δ t is the discrete time step, Ω( f i ) is the collision operator and F i is the implemented external forces term. Then, the collision operator in indicial form is as follows: where f eq i is the equilibrium function which is expressed by: The factors w i are given as: 4 9 , 1 9 , 0, 0, 1 9 , 1 36 , 0, 0, 1 36 , τ is the relaxation time without dimension. The nine discrete velocities are defined as follows: The collision term is expressed in the Multiple Relaxation Time model (MRT) D'Humières [37], where better stability is observed with a wide range of Prandtl Number values [38].
The flow field formulation becomes: where M is the projection matrix of f i and f eq i into the moment space. So, the expression of m = M f and m eq = M f eq are given by: The fluid density ρ, components moment j x and j y are the conserved quantities. The six other moments are non-conserved ones and are relaxed linearly in time namely: energy e, energy squared φ, energy flux in the two directions q x , q Y and diagonal/off-diagonal component of the strain-rate tensor p xx , p xy . The collision operator is carried out in the moment space and in indicial form: For the thermal model, the two-dimensional D2Q5 model with five velocities is used in this work. This model was chosen and validated by several authors in the literature [30] Fluids 2021, 6, 87 4 of 15 due to its simplicity and accuracy. The Boltzmann equation with multi-relaxation time can be written as: where g i (x, t) is distribution function of temperature, N is projection matrix of g i and g eq i into the moment space, and with the same procedure of temperature (first population) n = Ng. The transformation matrix N is given by: The boundary conditions, according to their macroscopic mathematical formulation are illustrated in Figure 1 for velocities: Before presenting the different results obtained, it is necessary to test the validity of our digital code. Thus, we compared the results of our simulations with those of the various studies carried out on the analysis of Rayleigh-Bénard convection in square cavities filled with air. Different discretization methods have been adopted by these reference studies. Table 1 quantitatively summarizes the values of the average Nusselt number obtained as a function of thermal gradient intensity imposed and characterized by the Rayleigh number. We note that our results show good agreement with those of the literature. The code is also successfully validated for the case of time-periodic temperature condition. Figure 2 shows a good agreement between the present hot and cold Nusselt numbers as function of time and those of Wang et al. [39].
The instantaneous and average Nusselt number Nu(t) at hot and cold wall (Y = 0 and Y = 1) is obtained as: and Nu c avg (t) = 1 0 Nu(X, t)dX (17) Before presenting the different results obtained, it is necessary to test the validity of our digital code. Thus, we compared the results of our simulations with those of the various studies carried out on the analysis of Rayleigh-Bénard convection in square cavities filled with air. Different discretization methods have been adopted by these reference studies. Table 1 quantitatively summarizes the values of the average Nusselt number obtained as a function of thermal gradient intensity imposed and characterized by the Rayleigh number. We note that our results show good agreement with those of the literature. The code is also successfully validated for the case of time-periodic temperature condition. Figure 2 shows a good agreement between the present hot and cold Nusselt numbers as function of time and those of Wang et al. [39]. (a) Comparison of the average Nusselt number at the hot wall with that of Wang et al. [39].
(b) Comparison of the average Nusselt number at the cold wall with that of Wang et al. [39]. In what follows, we present the influence of the various control parameters governing the natural convection problem. The different results are represented in terms of streamlines, isotherms and heat transfer rate over a time period. Finally, we analyze the periodicity of the convective regime, by means of phase portraits obtained by the heat exchange coefficients. In what follows, we present the influence of the various control parameters governing the natural convection problem. The different results are represented in terms of streamlines, isotherms and heat transfer rate over a time period. Finally, we analyze the periodicity of the convective regime, by means of phase portraits obtained by the heat exchange coefficients.

Results
The range of the Rayleigh number values is taken 10 3 ≤ Ra ≤ 10 6 . The amplitude of heating is taken between 0.0 and 0.8, and the working fluid is assumed to be air with Pr = 0.71.
In order to better study the behavior of the flow, we varied the amplitude for the three cases, 0.2, 0.5 and 0.8. Figure 3 represents the stream functions for the amplitude 0.2 as a function of the Rayleigh number for the four quarter-period f = 1/2πτ p = 1/3 × 10 −4 . Note that the stream function is obtained from the velocity integral and represents the fluid flow rate. When the Rayleigh number is equal to Ra = 10 4 , an invariant cell is observed for the four quarter-period.  This means that the modulation of cooling has little influence on the dynamic fields. The only difference is that the flow increases characterized by ψmax is increasingly larger than the general temperature gradient is maximum (in the heating period). Same case for the Rayleigh Ra = 10 5 with the presence of small vortices on the upper left and lower right side. As the thermal draft increases, for a Rayleigh number 10 6 , the flow loses its symmetry and subsequently becomes multicellular. For this Rayleigh value and a period of τp/2, the upper secondary cell grew further and took up most of the space of the cavity, in turn becoming counter-rotating bicellular. This spatial multi-cellular (mainly bicellular) competition is persistent throughout the period. Figures 4 and 5 represent the isocurrents for amplitudes 0.5 and 0.8 respectively, as This means that the modulation of cooling has little influence on the dynamic fields. The only difference is that the flow increases characterized by ψ max is increasingly larger than the general temperature gradient is maximum (in the heating period). Same case for the Rayleigh Ra = 10 5 with the presence of small vortices on the upper left and lower right side. As the thermal draft increases, for a Rayleigh number 10 6 , the flow loses its symmetry and subsequently becomes multicellular. For this Rayleigh value and a period of τp/2, the upper secondary cell grew further and took up most of the space of the cavity, in turn becoming counter-rotating bicellular. This spatial multi-cellular (mainly bicellular) competition is persistent throughout the period. Figures 4 and 5 represent the isocurrents for amplitudes 0.5 and 0.8 respectively, as a function of Rayleigh for the four quarter-period (with f = 1/3 × 10 −4 ). The flow structure for the Rayleigh 10 4 is essentially single-cell, with vortices on the upper right and lower left sides being very small. These findings are similar to the same Rayleigh with a smaller temperature modulation (Amp = 0.2). On the other hand, for Rayleigh equal to 10 5 , it is always the same type of flow (i.e., mainly single-cell) but the vortices of the corners are strongly present and influence the size of the main cell. The upper secondary cell still grows (Ra = 10 6 ) and the flow becomes predominantly counter-rotating bicellular and then becomes mainly monocellular again at the end of the period due to the decrease of the upper cell.  Finally, for the temporal step τp, this main cell will be oriented on the right side with one other cell on the lower left side. When the amplitude increases (i.e., Amp = 0.8, Figure  5), an essentially bicellular spatial competition is observed even for Ra = 10 5 and this competition persists even for the maximum Rayleigh of 10 6 .
The isotherms corresponding to the different Rayleigh number and the quarter-periods for Amp = 0.2 ( Figure 6) show that the heat distribution is consistent with the circulation of the fluid revealed by the stream functions. We also note that the isothermal lines are transported by the movement of fluid. For the Rayleigh 10 3 (not presented in this figure), the isotherms are stratified for all quarter-periods. The isotherms' distortion begins around Ra = 10 4 , evolving the form of a vortex. We notice that for Ra = 10 5 , isotherms are  When the Rayleigh reaches to 10 6 , we can clearly see the conformity of the temperature distribution with respect to the circulation of the fluid, and the distortion is always present for the four quarters of a period. It can be seen, for the last quarter of a period τp/2, that the temperature lines are concentrated much more on the cold wall compared to that of the opposite side; this is due to a cooling which is more intense. Finally, for the temporal step τp, this main cell will be oriented on the right side with one other cell on the lower left side. When the amplitude increases (i.e., Amp = 0.8, Figure 5), an essentially bicellular spatial competition is observed even for Ra = 10 5 and this competition persists even for the maximum Rayleigh of 10 6 .
The isotherms corresponding to the different Rayleigh number and the quarter-periods for Amp = 0.2 ( Figure 6) show that the heat distribution is consistent with the circulation of the fluid revealed by the stream functions. We also note that the isothermal lines are transported by the movement of fluid. For the Rayleigh 10 3 (not presented in this figure), the isotherms are stratified for all quarter-periods. The isotherms' distortion begins around Ra = 10 4 , evolving the form of a vortex. We notice that for Ra = 10 5 , isotherms are concentrated on the two horizontal walls for the four beats of the period. When the Rayleigh number reaches to 10 6 , we clearly see the conformity of the temperature distribution with respect to the fluid circulation, and distortion is always present for the four-quarter period. Figure 7 provides evolution of isotherms as function of Rayleigh number and the quarter-period for Amp = 0.8. Visibly, this figure shows that the heat distribution is consistent with the circulation of the fluid revealed by the stream functions. Additionally, we note that the isothermal lines are transported by the fluid flow. For Ra = 10 3 (not shown in the figure), isotherms are stratified for all quarter periods, and distortion begins around the Rayleigh value 10 4 in appearance of a vortex. We notice that for Raleigh equal to 10 5 , When the Rayleigh reaches to 10 6 , we can clearly see the conformity of the temperature distribution with respect to the circulation of the fluid, and the distortion is always present for the four quarters of a period. It can be seen, for the last quarter of a period τp/2, that the temperature lines are concentrated much more on the cold wall compared to that of the opposite side; this is due to a cooling which is more intense. When the Rayleigh reaches to 10 6 , we can clearly see the conformity of the temperature distribution with respect to the circulation of the fluid, and the distortion is always present for the four quarters of a period. It can be seen, for the last quarter of a period τp/2, that the temperature lines are concentrated much more on the cold wall compared to that of the opposite side; this is due to a cooling which is more intense. Figure 8 gives the temporal evolution of the Nusselt number at the cold wall for the last considered periods. The periodic oscillatory regime is established for most Rayleigh numbers. For high Rayleigh numbers, the cold Nusselt curves exhibit quasi-periodic shapes when the cooling modulation amplitude is small (Amp ≤ 0.5) and aperiodic shapes for high amplitudes (Amp > 0.5). We note that the more the amplitude increases, the more the heat transfer module intensifies. In the case of Rayleigh 10 6 , the curve loses its periodicity due to the appearance of unsteady phenomena. Note that when the amplitude increases, instabilities appear even for the intermediate Rayleigh numbers and lower than Ra < 10 6 . Additionally, it is apparent that for low modulation amplitudes, the curves have a quasi-symmetrical shape, whereas when the modulation increases, the physical behavior of the cavity is different between heating and cooling. It should also be noted that it is no longer the increase in the Rayleigh number that governs the growth of transfer rate.          Figure 10 illustrates the phase portraits of the normalized Nu c and Nu h (cold and hot Nusselt respectively) for different Ra, and for Amp = 0.4 and 0.8, respectively. It is well known that when boundary conditions are not modulated, hot and cold Nusselt numbers must be the same to ensure energy balance. In the case where the boundary conditions are modulated, the heat transfers on the hot and cold sides follow different evolutions [41][42][43][44]. The main observation is that for a Rayleigh of 10 3 , 10 4 and 10 5 , the limit cycle indicates that the regime is periodic; on the other hand for the Rayleigh 10 6 , the limit cycle is replaced by a cross cycle indicating the birth of natural instabilities and their addition to the pulsation imposed by the boundary condition. For Amp = 0.8, the limit and cross cycles are stretched diagonally showing that reaching the maximum value of Nu c induces the fall of Nu h to its minimum value and vice versa.

Conclusions
In this study, the Lattice Boltzmann method has been used in order to investiga Rayleigh-Bénard convection in square cavity submitted to a time-periodic cooling. T Rayleigh number value considered is between 10 3 and 10 6 , while the amplitude varies 0 to 0.8 and Prandtl value kept constant at 0.71. The flow state as well as the thermal fiel depends on the values of control parameters (Ra and Amp). For the dynamic field, t obtained results show that the flow structure changes from predominantly monocellul to predominantly counter-rotating bicellular flow for Rayleigh number values Ra= 10 6 an low values of the heating amplitude. This phenomenon was obtained for lower Rayleig

Conclusions
In this study, the Lattice Boltzmann method has been used in order to investigate Rayleigh-Bénard convection in square cavity submitted to a time-periodic cooling. The Rayleigh number value considered is between 10 3 and 10 6 , while the amplitude varies 0.2 to 0.8 and Prandtl value kept constant at 0.71. The flow state as well as the thermal fields depends on the values of control parameters (Ra and Amp). For the dynamic field, the obtained results show that the flow structure changes from predominantly monocellular to predominantly counter-rotating bicellular flow for Rayleigh number values Ra= 10 6 and low values of the heating amplitude. This phenomenon was obtained for lower Rayleigh values (Ra = 10 4 ) by increasing the value of the heating amplitude.
The analysis of the heat transfer show that for small values of amplitude heating, the averaged Nusselt curves follow a periodic evolution around an average corresponding to the formulation according to a constant cold temperature. An unsteady evolution is observed when the thermal draft increases; this unsteadiness has appeared for low Rayleigh numbers by increasing the value of the heating amplitude.

Conflicts of Interest:
The authors declare that they have no conflict of interest.