Analysis of unsteady mixed convection of Cu–water nanofluid in an oscillatory, lid-driven enclosure using lattice Boltzmann method

The unsteady physics of laminar mixed convection in a lid-driven enclosure filled with Cu–water nanofluid is numerically investigated. The top wall moves with constant velocity or with a temporally sinusoidal function, while the other walls are fixed. The horizontal top and bottom walls are, respectively, held at the low and high temperatures, and the vertical walls are assumed to be adiabatic. The governing equations along with the boundary conditions are solved through D2Q9 fluid flow and D2Q5 thermal lattice Boltzmann network. The effects of Richardson number and volume fractions of nanoparticles on the fluid flow and heat transfer are investigated. For the first time in the literature, the current study considers the mechanical power required for moving the top wall of the enclosure under various conditions. This reveals that the power demand increases if the enclosure is filled with a nanofluid in comparison with that with a pure fluid. Keeping a constant heat transfer rate, the required power diminishes by implementing a temporally sinusoidal velocity on the top wall rather than a constant velocity. Reducing frequency of the wall oscillation leads to heat transfer enhancement. Similarly, dropping Richardson number and raising the volume fraction of the nanoparticles enhance the heat transfer rate. Through these analyses, the present study provides a physical insight into the less investigated problem of unsteady mixed convection in enclosures with oscillatory walls.


Introduction
The fluid flow and heat transfer in the enclosures have been already studied extensively [1][2][3]. In the meantime, heat transfer in enclosures with moving boundaries has received considerable attention in recent years because of its growing application in the industry. Electronics cooling, solar collectors manufacturing, furnace engineering, lubrication technology and chemical equipment are some practical examples in which the physics of lid-driven enclosure is of importance [4].
Developments in nanotechnology applications have led to synthesis of nanoparticles, featuring high thermal conductivity [5,6]. Nanofluids, formed by uniformly dispersing nanoparticles through a pure base fluid, could significantly improve the heat transfer process in industrial equipment [7,8] and therefore attracted a great deal of interest; see for example [9][10][11]. Applying the advantages of nanoparticles to intensification of forced convection, some limited studies have been carried out on lid-driven enclosures filled with a nanofluid. These are presented in the following.
Karimipour et al. [12] investigated laminar mixed convection of copper-water nanofluid in a shallow inclined liddriven cavity. The upper hot wall moved with a constant velocity, while the other walls were maintained fixed and adiabatic. The governing equations were solved using lattice Boltzmann method (LBM). It was shown that higher Nusselt number could be achieved at higher inclined angles and volume fractions of nanoparticles mixed in the base fluid flow. This study demonstrated the excellent capability of LBM to simulate the problems including mixed convection of a nanofluid. LBM simulation of forced convection of water-Ag nanofluid in a square lid-driven cavity was analyzed by Dahani et al. [13]. The upper wall was pushed horizontally, while the other walls were fixed. The vertical walls held at a low temperature, while the horizontal walls were set adiabatic. A heated plate was fixed horizontally in the domain. The middle position of the heated plate was found to be the best location for optimizing the heat transfer.
A double lid-driven cavity containing a nanofluid of Cu-water was conducted using LBM by Rahmati et al. [14]. They considered mixed convection in the cavity with vertical wall, thermally sinusoidal heated. It was observed that the variation of the Richardson number and Nusselt number was in the opposite directions. Sheikholeslami et al. [15] studied nanofluid (alumina-water) convective heat transfer in a liddriven cubic with hot obstacles in the presence of a magnetic field using LBM. These authors considered the effects of Darcy and Reynolds numbers on the thermal boundary layer.
Zhou et al. [16] tackled the problem of double lid-driven enclosure containing nanofluid of alumina-water through using 3D LBM. The effects of nanoparticle volume fraction, Reynolds number and Richardson number on mixed convective heat transfer were taken into account. The best heat transfer characteristics were found for the case of downward moving of cold wall and upward moving for hot wall. It was shown that different heat source distributions, such as constant and sinusoidal function, significantly influence the temperature field and Nusselt number.
Through applying LBM, Sheikholeslami et al. [17] numerically conducted convection of heat and mass transfer in a permeable, non-Darcy cubic. The lower horizontal wall was moved, and a magnetic field was imposed upon the nanofluid of alumina-water. The domain contained a hot elliptic obstacle. Permeability, Reynolds number and Lorentz force were reported to dominate the physics of the problem. Further, multi-relaxation time LBM was applied to simulate a three-dimensional enclosure with lid-driven walls by Ghasemi and Siavashi [18]. Nusselt number was observed to increase at higher volume fraction values or Reynolds number.
Karimipour et al. [19] investigated periodic combined convection in a lid-driven cavity filled with a nanofluid using the finite volume method and SIMPLE algorithm to couple the pressure and velocity field. These authors reported smaller oscillations of mean Nusselt number at higher values of Richardson number. This study is apparently the only work on the lid-driven cavity with sinusoidal wall velocity in the literature.
The above review clearly depicts the shortage of studies on the problem of lid-driven enclosures filled with nanofluids by means of LBM. Given the importance of such configuration in industry, more investigations are required for understanding and optimization of the problem. Further, novel numerical methods, like LBM, are gradually spreading throughout the general field of computational fluid dynamics [20,21]. It follows that they should be used for classical problems, such as mixed convection in enclosures. Therefore, the current study aims to numerically investigate the mixed convection in a lid-driven enclosure filled with a nanofluid through using LBM. The rare boundary condition of sinusoidal motion of the wall is considered as a novel method of improving the thermohydraulic performance. Further, the power demand for the wall motion in periodic and constant velocity cases is compared to develop a comprehensive view of the problem.

Problem configuration and boundary conditions
The schematic geometry for the current problem is shown in Fig. 1. The height and width of the enclosure are represented, respectively, by H and L . The aspect ratio of the enclosure AR = L∕H is 3. Only the top wall moves, and the other walls are completely stationary. The top and bottom walls of the enclosure are correspondingly maintained at low ( T c ) and high temperatures ( T h ), while the vertical walls are assumed to be adiabatic. Gravitational acceleration is in the negative direction of y axis (Fig. 1).
Two distinct velocities for the top wall of the enclosure are assumed. These are a constant value, presented by u = u 0 , or a periodic sine function, shown as u = u 0 sin ( t) where u 0 is the amplitude of oscillations, t is time and denotes the angular frequency of motion. The physical problem and the corresponding numerical solution are therefore unsteady. This type of forcing has been already used to study nonlinear phenomena in many other fluid mechanical systems, such as low-density jets [22,23] and thermoacoustic systems [24,25].
The working fluid in the enclosure is water-Cu nanofluid, which is assumed Newtonian and incompressible. The flow is further considered laminar and single phase. The latter implies that the velocity of the fluid and solid particles is the same without any slip velocity. Boussinesq approximation is used to calculate the fluid density as a function of temperature [26,27].

Governing equations
The lattice Boltzmann method has been derived from the statistical view of the transport phenomena in the physical domain and solves the governing equations without any discretization [28]. Introducing a structured pattern for particles moving through the computation domain, the governing equations are determined. LBM comprises some relevant steps, namely streaming and collision [29,30], implemented in the current computational code, to solve the governing equations along with the boundary conditions. Once the convergence criterion is satisfied, the iterative process is terminated.  Fig. 2 The schematic of a D2Q9 and b D2Q5 networks As Fig. 2 illustrates, D2Q9 [31] and D2Q5 [32] networks are applied in the current work for the calculation of velocity and temperature, respectively.
At each time step, the density and velocity values for the D2Q9 network and the temperature upon D2Q5 network can be obtained from the following equations: in which c i is determined by [28] In Eq. (4), c is the base speed of the lattice Boltzmann, which is equal to unity [33]. The value of the equilibrium distribution function for the flow field and temperature are also obtained from [33] where V is the macroscopic velocity and W i is the weighting factor. On the flow field domain, W i takes the value of 4∕9 for i = 1 , 1∕9 for i = 2, 3, 4, 5 and 1∕36 at the other i values. However, the thermal field specifies the W i value of 1∕3 for i = 1 and 1∕6 when 2 ≤ i ≤ 5 . Since this work deals with mixed convection, the volume force of free convection should be defined in y-axis using the following relation [34]: In Eq. (6), c iy is the microscopic velocity of the particle in the y direction, g is the gravitational acceleration, is the volume expansion coefficient and T 0 is called the reference temperature, which equals the cold temperature ( T c ).

Non-dimensionalizing of the variables
Aiding the reference parameters, the governing hydrodynamic and thermal variables are turned into those of nondimensional form, by the following definitions. Next, the dimensionless parameters for more insightful physical analysis are stated in Table 1.
This process generalizes the results through reducing the wide variety of parameters and making the discussion more perceptive.
One of the most important non-dimensional variables is the Nusselt number. Therefore, spatially averaged Nusselt number is defined by where L is the length of the non-adiabatic horizontal walls. Because the current problem is time-dependent, the temporal averaging should be also performed to have a fully averaged Nusselt number in this problem. Considering as a temporal period of moving lid wall, the average Nusselt number becomes It should be pointed out that the Nusselt number is calculated after passing 50 temporal periods, due to wiping out of the probable irregular unsteadiness.

Mathematical description of boundary conditions
The boundary conditions are now presented in their nondimensional form. Dirichlet boundary conditions [35] are employed on the walls for velocity. On the top wall, the nonzero lid velocity and in the other walls zero velocity of no-slip condition are applied. These are expressed as the following: • The boundary conditions should be treated for the LBM. All the hydrodynamic distributed functions should be, therefore, specified on the boundaries. However, the unknown ones are those toward the domain, required to calculate. For example, f 6 , f 2 and f 9 need to be determined on the west wall of D2Q9 network, as shown in Fig. 3. This figure illustrates the schematic view of the problem with the direction of distributed functions on the walls. On each wall, the equations of density, x and y momentum, should be considered, as presented in Eqs. (18) to (20): There are three distribution functions as well as density on every wall, which are undetermined, and therefore, another equation is required in addition to Eqs. (18) to (20). The system of finding boundary conditions is closed using the non-equilibrium part of particle distribution normal to the boundary [36]. In this case, Eqs. (21) and (22) are, respectively, appropriate for the vertical and horizontal walls: It should be pointed out that the dimensionless velocity of the boundaries, i.e., u * and v * , is zero at all the walls except for the lid wall where it is determined by the constant or periodic functions, as stated earlier by Eq. (12).
Similar to what performed regarding hydrodynamic boundary conditions, unknown thermal distribution functions on the walls should be determined. However, as D2Q5 network is used for the energy equation, only a single unknown function, which is toward the domain, is required to be calculated. Aiding Eq. (3) and using the known temperature on the horizontal walls, the boundary condition equation is satisfied. On the vertical walls, all the distribution functions are assumed to be equal [33].

Nanofluid properties
As stated earlier, single phase assumption of nanofluid mixture is considered in the current study. This approach uses some mathematical relations to determine the thermophysical properties of the nanofluid. Nanofluid density can be determined by the following equation [37]: where subscribes f, s and nf denote fluid, solid and nanofluid, respectively. and are the volume fraction of the nanoparticles and density, correspondingly. The heat capacity, c p , and the volume expansion coefficient, , for the nanofluid are also defined as follows [37]: The nanofluid dynamic viscosity, nf , is calculated using the Brinkman model [38]: The effective conductivity coefficient of the nanofluid is defined by the Maxwell model [39], which is The thermo-physical properties of the current nanofluid component (copper nanoparticles and pure water) at the reference temperature of 298 K are presented in Table 2.

Grid independency
Grid independency test was performed to choose the optimal mesh size. The mesh used in the current study is structural with a uniform size in each direction. The parameters value applied during the grid independency test is presented in Table 3. Further, the non-dimensional time step, Δt * , was assumed to be unity, for restricting the CFL 1 number below those recommended for the similar problems [40].
Four mesh sizes of 40 × 120, 50 × 150, 60 × 180 and 70 × 210 on the x-y domain were examined. The nondimensional temperature versus time at the midpoint of the   enclosure, located at x * = y * = 0.5 , is presented in Fig. 4 for various mesh sizes. This figure shows that the mesh size of 60 × 180 is more precise and economic than the others and therefore was chosen for the rest of simulations.

Validation of the computational code
In order to validate the in-house computational code, the current results were compared to the other published results. The current results for dimensionless velocity and the local Nusselt number on the mid-vertical centerline of a two-sided lid-driven cavity filled with Cu-water nanofluid are compared with those of Ref. [38], respectively, as shown in Fig. 5a, b. Clearly, the two sets of data are in very good agreement. Finally, the current results are also compared with those of Karimipour et al. [19], as shown in Fig. 6. This figure includes the dimensionless temperature for the nondimensional time that matches the temporal period of lid wall motion at the vertical centerline. A high degree of agreement is observed in this figure. These comparisons can potentially prove that the current in-house code can completely capture the physics of fluid flow and heat transfer in the transient lid-driven enclosure.

Results and discussion
In this section, the results of mixed convection in a sinusoidal lid-driven enclosure filled with copper-water nanofluid are presented for Gr = 10 4 Pr = 6.2 and L∕H = 3 . The time-dependent motion of the top wall causes the temporal variation of the parameters, as it will be shown in this section. Further, the effects of changes in Richardson number and nanoparticles volume fraction and the lid wall period on the fluid flow and heat transfer are investigated. Figure 7 shows the influences of variations in Richardson number on the average Nusselt number. According to this figure, decreasing the numerical value of Richardson number leads to increases in the average Nusselt number, as the lid velocity rises. This figure shows that decrement of R from 10 to 0.1 makes the averaged Nu number increase by the factor of 1.8. At the lowest lid velocity ( Ri = 0.1 ), variation of the average Nu number demonstrates a decreasing trend with growing frequencies, representing a low-pass-filter trend in this problem. However, the average Nusselt number shows steady-like variation at the highest value of Ri, emanated from the lowest velocity of the lid wall.
The non-dimensional horizontal velocity at the vertical mid-centerline in one period of lid wall moving is illustrated for two values of Ri and zero nanoparticles volume fraction in Fig. 8. The top lid velocity is positive when the nondimensional time is lower than ∕2 . In this case, the shear stress is applied in the positive direction of x * . However, this trend is reversed for the non-dimensional times greater than ∕2 . The velocity profile of non-dimensional time that is equal to 0 and ∕2 and further ∕4 and 3 ∕4 is mirror image with respect to each other as the simulation is performed under a stationary unsteady condition. The variation of the velocity is negligible around the core of the enclosure compared to that in the vicinity of the wall regions. Farther away from the top wall, the lid velocity penetration becomes less significant. Yet, as the non-dimensional frequency of the lid oscillation decreases, the velocity effects of the lid penetrate deeper into the fluid. Lower velocity at the vicinity of the bottom wall zone renders a weaker shear stress at that region. Figure 9, drawn for the same conditions as in Fig. 8, shows the non-dimensional temperature inside the enclosure. The temperature profile is varied between the minimum and maximum values, which is 1 and − 1. The symmetric profile at Ri = 10 shows the uniform free convection cell through the enclosure. Further, higher Richardson number makes the buoyancy effects stronger and this is the reason for the extension of the higher temperature zone close to the horizontal walls region. The temperature profile does not change by time passing after the beginning ( t * = 0 ). This shows a fast thermal response of the system, leading to stationary conditions for the temperature just after the lid motion. Figures 10-13 show the streamline and isotherms for one period of lid wall oscillation and for various Richardson numbers and volume fractions. At the even quarters of temporal period and Ri = 10 , most of the enclosure is filled with a large cell for zero concentration of nanoparticles. Further, one or two small cells may be found in the domain. However, this pattern changes by decreasing the Richardson number or increasing Reynolds number. In this case, the number of recirculating cells is increased. The different cells arrangement between the cases of = 0% and = 4% arises from the difference between dynamic viscosity. At odd quarters of time period, three recirculating cells occupy the enclosure at Ri = 10 , while one big wake as well as another small cell covers the domain at Ri = 0.1 . Isotherms show the conduction of heat on the non-adiabatic walls, represented by the dense and parallel lines. This type of heat transfer tends to magnify at higher Richardson numbers. The thickness of thermal boundary layer is found to be lower at lower Ri , as a result of intensification of forced convection by increasing the lid wall velocity. This also causes thickness reduction in the hydrodynamic boundary layer. The patterns of the isotherms and streamlines are similar at t * = 0 , t * = ∕2 and t * = . Further, the distribution of the isotherms and streamlines of t * = 3 ∕4 is the same as those of t * = ∕4 , rendering a stationary temporal behavior.
Variations of Nusselt number versus three periods of lid wall motion at Ri = 0.1 and 10 and for various nanoparticles volume fractions are illustrated in Fig. 14. The Nusselt number is periodically changed due to the sinusoidal velocity of the lid wall. This clearly shows that the heat transfer behavior of the enclosure is time-dependent and varies according to the frequency of the wall moving. Interestingly, about two periods of Nusselt number variation is set on one period of lid wall motion. This is a clear manifestation of the nonlinear behavior of the system in which the frequency of the input forcing (lid oscillation) is different to that of the output (Nusselt number). Reducing Ri from 10 to 0.1 causes the Nusselt number increase by a factor of two. This is physically conceivable as the forced convection is stronger than the free convection. Increasing the volume fraction makes an increment in convection heat transfer [41]. However, its growth rate is not restricted by increasing Ri, in contrast to the previous studies [38]. Therefore, as an important result, wall oscillations can compensate this weakness in heat transfer. Figure 15 depicts a comparison of the average Nusselt number for different values of Ri and the volume fraction of nanoparticles. Clearly, a reduction in Ri results in higher rates of heat transfer. This could be readily explained by noting that forced convection of heat is strengthened at lower values of Ri. Enhancement of Nusselt number is also achieved by increasing the concentration of nanoparticles, as it boosts the thermo-physical properties of the nanofluid. This figure shows that heat transfer is partly subsided by imposing a sinusoidal motion on the lid wall, especially at lower Richardson numbers. It is inferred from Fig. 15 that throughout the considered range of the parameters, the effects of Ri are stronger than that of the volume fraction of nanoparticles. The non-dimensional frequencies in Fig. 15 are those that lead to obtaining the maximum Nusselt number. Intensifying the non-dimensional frequency of the lid oscillation for the cases with higher velocity ( Ri ≤ 1 ) may be another reason for the decline of Nusselt number. This behavior is further expected from the lowpass filter response of the system.
The mechanical power required to move the lid wall is calculated for constant or oscillatory wall velocity in this study. Importantly, it is noted that despite its practical importance, so far, this aspect of the problem has not been The imposed force on the lid wall is computed using shear stress, w , as the following relation shows:  The average power on one oscillation period gives Figure 16 illustrates the average heat transfer rate out of the top wall against the average power demand to move the enclosure's wall. As expected, increasing heat transfer magnifies the power consumption to drive the wall. The difference between two cases of constant and oscillatory velocity of the lid wall is small at lower power demands, while it grows to a higher extent as power input rises. Figure 16 shows that for a constant heat transfer rate, the enclosure with constant lid velocity needs more input power than that of the periodic velocity case. This indicates that for a specific thermal gain, the enclosure of the sinusoidal wall velocity is more efficient. Lower power demand for the oscillatory wall may be induced by lower average velocity in the temporal period compared to the constant lid velocity. The variation of heat transfer rate through the graph tends to decline by increasing the power demand. This implies that as the power input of the system increases, the heat transfer system takes farther distance from its optimal operation.
The average rate of heat transfer versus power demand for two values of volume fraction for the case of periodic wall velocity is depicted in Fig. 17. The power demand for a constant heat transfer rate grows by increasing the volume fraction. This emanates from the fact that a nanofluid flow involves higher shear stress due to more dynamic viscosity in comparison with the base fluid flow. In keeping with the findings of other thermohydraulic studies of nanofluid flows [42][43][44][45][46] and those on other types of fluids [47][48][49][50], the power demand grows for the cases of higher heat transfer.

Conclusions
Using the lattice Boltzmann method, mixed heat convection in a two-dimensional lid-driven enclosure filled with nanofluid was simulated. The lid velocity was assumed to be oscillatory and sinusoidal. A range of Richardson number ( 0.1 to 10 ) and volume fraction of nanoparticle ( 0 to 4% ) was considered. For the first time in the literature, isotherms, streamlines and Nusselt number as well as the power demand to move the lid wall were studied. The main findings of this study are summarized as follows.
• The enclosure with sinusoidal lid wall velocity was demonstrated to act as a classical low-pass filter, which responds most strongly to low-frequency oscillations and filters out higher frequencies. This was the reason for decline of Nusselt number as the non-dimensional frequency increased. • Assuming a constant heat transfer rate from the enclosure, the case with oscillatory lid wall velocity showed a lower demand of input power in comparison with that with constant lid wall speed. • The power required to oscillate the lid wall grew by increasing the concentration of nanoparticles, due to increases in the fluid viscosity. • Keeping Grashof number at a constant value, the Nusselt number increased with decreasing Richardson number, due to magnification of forced convection.
Increasing the volume fraction of nanoparticles also appeared to boost the heat transfer rate.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.