*Article* **Time-Periodic Cooling of Rayleigh–Bénard Convection**

#### **Lyes Nasseri <sup>1</sup> , Nabil Himrane <sup>2</sup> , Djamel Eddine Ameziani <sup>1</sup> , Abderrahmane Bourada <sup>3</sup> and Rachid Bennacer 4,\***


**Abstract:** 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<sup>3</sup> <sup>≤</sup> Ra <sup>≤</sup> <sup>10</sup><sup>6</sup> 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.

**Keywords:** Rayleigh–Bénard convection; time periodical cooling; Lattice Boltzmann method

#### **1. 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–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–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.

**Citation:** Nasseri, L.; Himrane, N.; Ameziani, D.E.; Bourada, A.; Bennacer, R. Time-Periodic Cooling of Rayleigh–Bénard Convection. *Fluids* **2021**, *6*, 87. https://doi.org/ 10.3390/fluids6020087

Academic Editor: Marcello Lappa

Received: 19 January 2021 Accepted: 12 February 2021 Published: 16 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–15]. Other researchers were interested in the Rayleigh-Bénard convection in viscoplastic fluids leading to numerical and experimental investigations [14–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–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–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.

#### **2. 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]:

$$f\_i(\mathbf{x} + e\_l \delta\_{l\prime} t + \delta\_{l\prime}) - f\_i(\mathbf{x}, t) = \Omega(f\_i) + \delta\_l F\_i \tag{1}$$

where *f<sup>i</sup>* is the distribution function with velocity *e<sup>i</sup>* at lattice node *x* at time *t*, *δ<sup>t</sup>* is the discrete time step, *Ω*(*fi*) is the collision operator and *F<sup>i</sup>* is the implemented external forces term. Then, the collision operator in indicial form is as follows:

$$f\_l(\mathbf{x} + e\_l \delta\_l, t + \delta\_l) - f\_l(\mathbf{x}, t) = -1/\pi \left[ f\_l(\mathbf{x}, t) - f\_l^{eq}(\mathbf{x}, t) \right] + \delta\_l F\_l \tag{2}$$

where *f eq i* is the equilibrium function which is expressed by:

$$f\_i^{eq} = w\_i \rho \left[1 + 3e\_i v + 9(e\_i v)^2 / 2 - 3v^2 / 2\right] \tag{3}$$

The factors w<sup>i</sup> are given as: <sup>n</sup> 4 9 , 1 9 , 0, 0, <sup>1</sup> 9 , 1 <sup>36</sup> , 0, 0, <sup>1</sup> 36o , *τ* is the relaxation time without dimension. The nine discrete velocities are defined as follows:

$$e\_i = \begin{cases} (0,0); \; i = 0\\ \csc[\cos((i-1)\pi/2), \sin((i-1)\pi/2)]; \; i = 1,2,3,4\\ \sqrt{2c}[\cos((2i-9)\pi/2), \sin((2i-9)\pi/2)]; \; i = 5,6,7,8 \end{cases} \tag{4}$$

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].

$$
\Omega = -M^{-1}\mathbb{C}[m\_i(\mathbf{x}, t) - m\_i^{\epsilon q}(\mathbf{x}, t)] \tag{5}
$$

The flow field formulation becomes:

$$f\_i(\mathbf{x} + e\_l \delta\_l, t + \delta\_l) - f\_i(\mathbf{x}, t) = -M^{-1} \mathbb{C}[m\_i(\mathbf{x}, t) - m\_i^{eq}(\mathbf{x}, t)] + M^{-1} \delta\_l (1 - \mathbb{C}/2)D \tag{6}$$

where *M* is the projection matrix of *f<sup>i</sup>* and *f eq i* into the moment space. So, the expression of *m* = *M f* and *meq* = *M f eq* are given by:

$$
\begin{pmatrix}
\rho\\ \epsilon\\ \Phi\\ j\_{\mathbf{x}}- (\delta\_{l}/2)\rho F\_{\mathbf{x}}\\ q\_{\mathbf{x}}\\ j\_{\mathbf{y}}- (\delta\_{l}/2)\rho F\_{\mathbf{y}}\\ q\_{\mathbf{y}}\\ p\_{\mathbf{x}\mathbf{y}}\\ p\_{\mathbf{xy}}
\end{pmatrix} = 
\begin{bmatrix}
1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\ -4 & -1 & -1 & -1 & -1 & 2 & 2 & 2 & 2\\ 4 & -2 & -2 & -2 & -2 & 1 & 1 & 1 & 1\\ 0 & 1 & 0 & -1 & 0 & 1 & -1 & -1 & 1\\ 0 & -2 & 0 & 2 & 0 & 1 & -1 & -1 & 1\\ 0 & 0 & 1 & 0 & -1 & 1 & 1 & -1 & -1\\ 0 & 0 & -2 & 0 & 2 & 1 & 1 & -1 & -1\\ 0 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1
\end{bmatrix}
\begin{pmatrix}
f\_{0} \\ f\_{1} \\ f\_{2} \\ f\_{3} \\ f\_{4} \\ f\_{5} \\ f\_{6} \\ f\_{7} \\ f\_{8} \\ f\_{9}
\end{pmatrix}.
\tag{7}
$$

The fluid density *ρ*, components moment *j<sup>x</sup>* and *j<sup>y</sup>* 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 *qx*, *q<sup>Y</sup>* and diagonal/off-diagonal component of the strain-rate tensor *pxx*, *pxy*. The collision operator is carried out in the moment space and in indicial form:

$$m\_i^\*(\mathbf{x}, t) = m\_i(\mathbf{x}, t) - \mathbb{C}[m\_i(\mathbf{x}, t) - m\_i^{eq}(\mathbf{x}, t)] \tag{8}$$

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] due to its simplicity and accuracy. The Boltzmann equation with multi-relaxation time can be written as:

$$\mathcal{g}\_i(\mathbf{x} + e\_i \delta\_{\mathbf{t}\prime} t + \delta\_{\mathbf{t}}) - \mathcal{g}\_i(\mathbf{x}, t) = -\mathbf{N}^{-1} \mathcal{E} [n\_i(\mathbf{x}, t) - n\_i^{eq}(\mathbf{x}, t)] \tag{9}$$

where *gi*(*x*, *t*) is distribution function of temperature, *N* is projection matrix of *g<sup>i</sup>* 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:

$$
\begin{pmatrix} n\_0 \\ n\_1 \\ n\_2 \\ n\_3 \\ n\_4 \end{pmatrix} = \begin{bmatrix} 1 & 1 & 1 & 1 & 1 \\ 0 & 1 & 0 & -1 & 0 \\ 0 & 0 & 1 & 0 & -1 \\ -4 & 1 & 1 & 1 & 1 \\ 0 & 1 & -1 & 1 & -1 \end{bmatrix} \begin{pmatrix} \text{g} \\ \text{g} \\ \text{g} \\ \text{g}\_3 \\ \text{g}\_4 \end{pmatrix} \tag{10}
$$

The boundary conditions, according to their macroscopic mathematical formulation are illustrated in Figure 1 for velocities:

$$\mathcal{U}(0,Y,t) = \mathcal{U}(1,Y,t) = \mathcal{U}(X,0,t) = \mathcal{U}(X,1,t) = 0 \tag{11}$$

$$V(0, Y, t) = V(1, Y, t) = V(X, 0, t) = V(X, 1, t) = 0 \tag{12}$$

$$\mathbf{u}\_{\mathbf{c}} = \mathbf{A}\mathbf{m}\mathbf{p} \times \sin\left(2\pi f \mathbf{r}t\right),$$

**Figure 1.** Physical model. **Figure 1.** Physical model.

> Before presenting the different results obtained, it is necessary to test the validity of Temperature boundary conditions:

$$\theta(X, 1, t) = Amp \cdot \sin(2\pi ft) \tag{13}$$

$$\theta(\mathbf{X}, \mathbf{0}, t) = 1 \tag{14}$$

$$\left. \frac{\partial \theta(X, Y, t)}{\partial X} \right|\_{X=0} = 0 \text{ and } \left. \frac{\partial \theta(X, Y, t)}{\partial X} \right|\_{X=1} = 0 \tag{15}$$

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:

$$\mathrm{Nu}\_{\hbar}(X,0,t) = -\frac{\partial \theta(X,Y,t)}{\partial Y}\bigg|\_{Y=0} \quad \text{and} \, \mathrm{Nu}\_{\hbar\_{\text{avg}}}(t) = \int\_{0}^{1} \mathrm{Nu}(X,t)dX \tag{16}$$

**103 104 105 106**

Present Work 1.0035 2.1502 3.912 6.321 Ourtatani et al. [13] 1.0004 2.158 3.910 6.309 Turan [14] 1.0000 2.154 3.907 6.309 Bouabdallah et al. [40] 1.0000 2.2000 3.900 6.400

$$\mathrm{Nu}\_{\mathfrak{c}}(X, 1, t) = -\frac{\partial \theta(X, Y, t)}{\partial Y}\Big|\_{Y=1} \text{ and } \mathrm{Nu}\_{\mathfrak{c}\_{\mathfrak{w}\_{\mathfrak{X}}}}(t) = \int\_0^1 \mathrm{Nu}(X, t)dX \tag{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].

**Figure 2.** Time history of Nusselt numbers. Ra = 106, Pr = 6.2, τp = 0.1 and Amp = 0.8. **Figure 2.** Time history of Nusselt numbers. Ra = 10<sup>6</sup> , Pr = 6.2, τp = 0.1 and Amp = 0.8.

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 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

The range of the Rayleigh number values is taken 103 ≤ Ra ≤ 106. The amplitude of heating is taken between 0.0 and 0.8, and the working fluid is assumed to be air with

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 <sup>4</sup> 1 2 1 3 10

. 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 = 104, an invariant cell is ob-

πτ<sup>−</sup> = =⋅ *<sup>p</sup> f*

exchange coefficients.

served for the four quarter-period.

**3. Results** 

Pr 0.71 = .

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.

#### **3. Results**

The range of the Rayleigh number values is taken 10<sup>3</sup> <sup>≤</sup> Ra <sup>≤</sup> <sup>10</sup><sup>6</sup> . 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 *<sup>f</sup>* <sup>=</sup> 1/2*πτ<sup>p</sup>* <sup>=</sup> 1/3 <sup>×</sup> <sup>10</sup>−<sup>4</sup> . 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<sup>4</sup> , an invariant cell is observed for the four quarter-period. *Fluids* **2021**, *6*, 87 7 of 16

**Figure 3.** Evolution of streamlines over a period as a function of the Rayleigh number (Amp = 0.2 and |f| = 0.33 × 10<sup>−</sup>4). **Figure 3.** Evolution of streamlines over a period as a function of the Rayleigh number (Amp = 0.2 and |f| = 0.33 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ).

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 = 105 with the presence of small vortices on the upper left and lower right side. As the thermal draft increases, for a Rayleigh number 106, the flow loses its symmetry 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<sup>5</sup> with the presence of small vortices on the upper left and lower right side. As the thermal draft increases, for a Rayleigh number 10<sup>6</sup> , the flow loses its

and subsequently becomes multicellular. For this Rayleigh value and a period of τp/2, the

becoming counter-rotating bicellular. This spatial multi-cellular (mainly bicellular) com-

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 104 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 105, 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 = 106) 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

petition is persistent throughout the period.

upper cell.

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 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ). The flow structure for the Rayleigh 10<sup>4</sup> 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<sup>5</sup> , 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<sup>6</sup> ) 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. *Fluids* **2021**, *6*, 87 8 of 16

**Figure 4.** Evolution of streamlines over a period as a function of the Rayleigh number (Amp = 0.5 and |f| = 0.33 × 10<sup>−</sup>4). **Figure 4.** Evolution of streamlines over a period as a function of the Rayleigh number (Amp = 0.5 and |f| = 0.33 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ).

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 = 105 and this com-

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 circula-

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 = 103(not shown in the figure), isotherms are stratified for all quarter periods, and distortion begins around the Rayleigh value 104 in appearance of a vortex. We notice that for Raleigh equal to 105,

are transported by the movement of fluid. For the Rayleigh 103 (not presented in this figure), the isotherms are stratified for all quarter-periods. The isotherms' distortion begins around Ra = 104, evolving the form of a vortex. We notice that for Ra = 105, isotherms are concentrated on the two horizontal walls for the four beats of the period. When the Rayleigh number reaches to 106, 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.

petition persists even for the maximum Rayleigh of 106.

10<sup>−</sup>4).

Ra = 104

min max θ θ

times of the period.

**Figure 5.** Evolution of streamlines "Ψ" over a period as a function of the Rayleigh number (Amp = 0.8 and |f| = 0 × 33 **Figure 5.** Evolution of streamlines "Ψ" over a period as a function of the Rayleigh number (Amp = 0.8 and |f| = 0 <sup>×</sup> 33 10−<sup>4</sup> ).

When the Rayleigh reaches to 106, 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 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<sup>5</sup> and this competition persists even for the maximum Rayleigh of 10<sup>6</sup> .

the temperature lines are found concentrated on the two horizontal walls for the four

τ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. 0 or 4/4τp 1/4τp 2/4τp 3/4τ<sup>p</sup> 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<sup>3</sup> (not presented in this figure), the isotherms are stratified for all quarter-periods. The isotherms' distortion begins around Ra = 10<sup>4</sup> , evolving the form of a vortex. We notice that for Ra = 10<sup>5</sup> , isotherms are concentrated on the two horizontal walls for the four beats of the period. When the Rayleigh number reaches to 10<sup>6</sup> , 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.

 / = 0.010 / 0.996 min max θ θ/ =−0.190 / 0.996 min max θ θ/ =−0.017 / 0.996 min max θ θ/ = 0.158 / 0.996 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<sup>3</sup> (not shown in the figure), isotherms are stratified for all quarter periods, and distortion begins around the Rayleigh value 10<sup>4</sup> in appearance of a vortex. We notice that for Raleigh equal to 10<sup>5</sup> ,

Ra = 104

Ra = 105

Ra = 106

min max θ θ

min max θ θ

min max θ θ

θ θ

/ =0.430/ 0.992 min max

θ θ

/ =−0.333/0.993 min max

Ra = 104

Ra = 105

Ra = 106

*Fluids* **2021**, *6*, 87 9 of 16

0 or 4/4τp 1/4τp 2/4τp 3/4τ<sup>p</sup>

min max Ψ Ψ =− / 0.001/ 6.905 min max Ψ Ψ =− / 0.001/ 5.158 min max Ψ Ψ =− / 0.001/ 6.223 min max Ψ Ψ =− / 0.001/ 8.179

min max Ψ Ψ =− / 1.854 /14.79 min max Ψ Ψ =− / 6.431/ 22.01 min max Ψ Ψ =− / 2.941/ 24.03 min max Ψ Ψ =− / 0.270 / 24.74

min max Ψ Ψ =− / 13.87 / 54.88 min max Ψ Ψ =− / 60.12 / 40.88 min max Ψ Ψ =− / 53.42 / 57.00 min max Ψ Ψ =− / 85.28 / 6.686

to that of the opposite side; this is due to a cooling which is more intense.

**Figure 5.** Evolution of streamlines "Ψ" over a period as a function of the Rayleigh number (Amp = 0.8 and |f| = 0 × 33

times of the period.

the temperature lines are found concentrated on the two horizontal walls for the four

the temperature lines are found concentrated on the two horizontal walls for the four times of the period. 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

When the Rayleigh reaches to 106, we can clearly see the conformity of the temperature distribution with respect to the circulation of the fluid, and the distortion is always

**Figure 6.** Evolution of isotherms over a period as a function of the Rayleigh number (Amp = 0.2 and |f| = 0.33 × 10<sup>−</sup>4). **Figure 6.** Evolution of isotherms over a period as a function of the Rayleigh number (Amp = 0.2 and |f| = 0.33 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ).

 0 or 4/4τp 1/4τp 2/4τp 3/4τ<sup>p</sup> When the Rayleigh reaches to 10<sup>6</sup> , 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.

 / = 0.070 / 0.996 min max θ θ/ = 0.031/ 0.996 min max θ θ/ =−0.771/ 0.995 min max θ θ/ =−0.240 / 0.996 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<sup>6</sup> , 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<sup>6</sup> . 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.

/ =−0.779/0.996 min max

/ =−0.782 / 0.993 min max

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

θ θ

θ θ

/ =−0.163 / 0.995

/ =−0.146/ 0.995

/ =0.024/ 0.995 min max

**Figure 7.** Evolution of isotherms over a period as a function of the Rayleigh number (Amp = 0.8 and |f| = 0.33 × 10<sup>−</sup>4).

θ θ

θ θ

periods.

Nu

 Amp 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9

Ra = 105

Ra = 106

min max θ θ

min max θ θ

θ θ

/ = 0.008 / 0.995 min max

θ θ

/ = 0.011/ 0.995 min max

/ =−0.191/ 0.995 min max

/ =−0.192 / 0.994 min max

**Figure 6.** Evolution of isotherms over a period as a function of the Rayleigh number (Amp = 0.2 and |f| = 0.33 × 10<sup>−</sup>4).

θ θ

θ θ

/ = 0.000 / 0.994 min max

/ = 0.000 / 0.995 min max

θ θ

θ θ

/ = 0.201/ 0.994

/ = 0.200 / 0.994

**Figure 7.** Evolution of isotherms over a period as a function of the Rayleigh number (Amp = 0.8 and |f| = 0.33 × 10<sup>−</sup>4). **Figure 7.** Evolution of isotherms over a period as a function of the Rayleigh number (Amp = 0.8 and |f| = 0.33 <sup>×</sup> <sup>10</sup>−<sup>4</sup> ).

**Figure 8.** Evolution of the Nusselt number "Nuc" (at the top wall) as a function of time for different Ra over the last two **Figure 8.** *Cont*.

Nu

 Amp 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9

840000 850000 860000 870000 880000 890000

Ra = 10<sup>4</sup>

t

840000 850000 860000 870000 880000 890000

Ra = 10<sup>3</sup>

t

**Figure 8.** Evolution of the Nusselt number "Nuc" (at the top wall) as a function of time for different Ra over the last two periods. **Figure 8.** Evolution of the Nusselt number "Nuc" (at the top wall) as a function of time for different Ra over the last two periods. **Figure 8.** Evolution of the Nusselt number "Nuc" (at the top wall) as a function of time for different Ra over the last two periods.

 **Figure 9.** Evolution of the Nusselt number "Nuc" (at the top wall) as a function of time for different Amp over the last two periods. **Figure 9.** Evolution of the Nusselt number "Nuc" (at the top wall) as a function of time for different Amp over the last two periods.

t

Figure 10 illustrates the phase portraits of the normalized Nuc and Nuh (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– 44]. The main observation is that for a Rayleigh of 103, 104 and 105, the limit cycle indicates that the regime is periodic; on the other hand for the Rayleigh 106, the limit cycle is replaced by a cross cycle indicating the birth of natural instabilities and their addition to the

t

Figure 10 illustrates the phase portraits of the normalized Nu<sup>c</sup> and Nu<sup>h</sup> (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–44]. The main observation is that for a Rayleigh of 10<sup>3</sup> , 10<sup>4</sup> and 10<sup>5</sup> , the limit cycle indicates that the regime is periodic; on the other hand for the Rayleigh 10<sup>6</sup> , 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<sup>c</sup> induces the fall of Nu<sup>h</sup> to its minimum value and vice versa. *Fluids* **2021**, *6*, 87 13 of 16 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 Nuc induces the fall of Nuh to its minimum value and vice versa.

**Figure 10.** Phase plot for the normalized Nuc and Nuh for different Ra (Amp = 0.4 and 0.8). **Figure 10.** Phase plot for the normalized Nu<sup>c</sup> and Nu<sup>h</sup> for different Ra (Amp = 0.4 and 0.8).

#### **4. Conclusions 4. 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 103 and 106, 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 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<sup>3</sup> and 10<sup>6</sup> , 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

values (Ra = 104) by increasing the value of the heating amplitude.

obtained results show that the flow structure changes from predominantly monocellular to predominantly counter-rotating bicellular flow for Rayleigh number values Ra= 106 and low values of the heating amplitude. This phenomenon was obtained for lower Rayleigh

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

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<sup>6</sup> and low values of the heating amplitude. This phenomenon was obtained for lower Rayleigh values (Ra = 10<sup>4</sup> ) 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.

**Author Contributions:** Conceptualization, D.E.A.; Formal analysis, L.N., N.H., D.E.A. and R.B.; Investigation, L.N. and N.H.; Methodology, L.N. and N.H.; Supervision, R.B.; Writing—original draft, L.N.; Writing—review & editing, N.H., D.E.A. and A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** No data available.

**Acknowledgments:** The authors acknowledge the financial support of the General Direction of Scientific Research and Technological Development (DGRSDT), Algeria.

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

#### **Abbreviations**


#### **Subscript**


#### **Greek symbols**


#### **References**

