Numerical Study of the Effect of Viscous Heat Dissipation and Compression Work on Microscale Rayleigh – Bénard Convection Based on a Coupled Thermal Lattice Boltzmann Method

In the present work, an improved double-distribution-function thermal lattice Boltzmann method (LBM) is developed for analyzing the effect of viscous heat dissipation and compression work on microscale Rayleigh–Bénard convection. In the proposed method a temperature change is introduced into the LB momentum equation in the form of a momentum source to realize the coupling between the momentum and the energy fields; two sets of evolution equations are established, one for the mass and momentum conservation and the other for the total energy that incorporates viscous heat dissipation and compression work. Numerical results show that the effect of viscous heat dissipation and compression work on the temperature distribution, flow distribution, and average Nusselt number at some Rayleigh numbers and aspect ratios is significant.


INTRODUCTION
Rayleigh-Bénard convection [1] is a thermally induced flow between two horizontal plates that are heated from below.As one of benchmark problems of hydrodynamic instability, Rayleigh-Bénard convection is of considerable scientific and engineering importance.Moreover, Rayleigh-Bénard convection involves in various technological applications such as nuclear reactor insulation, solar collectors, crystal growth in liquids, heat exchangers, and cooling of electronic equipment, and particularly the cooling of electronic circuits in microelectromechanical systems [2].Therefore, much effort has been spent on investigating Rayleigh-Bénard convection.Zhou et al. studied the geometric and physical properties of thermal plumes in turbulent Rayleigh-Bénard convection [3].A high-resolution measurement of the velocity and temperature boundary layers in Rayleigh-Bénard convection was conducted by Sun et al. with applying the particle image velocimetry technique [4].Sharif and Mohammad investigated the natural convection within cavities with different inclination angles and cavity aspect ratios [5].Valencia et al. measured the velocity field and timeaveraged flow structures in a cubical cavity that was heated from below and cooled from above [6].Kao et al. investigated the nonlinear phenomena of two-dimensional (2D) natural convection in enclosed rectangular cavities and systematically examined the relationship ___________________________________ *Corresponding Author: sswei@sdu.edu.cn between the Nusselt number (Nu) and the reference Rayleigh number (Ra) [7].Chakraborty and Chatterjee studied the Rayleigh-Bénard convection occurring in directional solidification by using a novel hybrid lattice Boltzmann method (LBM) [8].
The aforementioned investigations showed that the formation of flow structures due to Rayleigh-Bénard convection are affected by the boundary conditions, aspect ratio of the cavity, and direction of the gravitational effect.However, few of the studies, except that presented in [9], have considered the effect of viscous heat dissipation and compression work on microscale Rayleigh-Bénard convection.In the model developed in Ref. 9 the temperature equilibrium distribution function always takes negative values, and this model is more complicated than other double-distribution-function (DDF) models, even when the viscous heat dissipation and compression work are negligible.
With a decrease of the system size, the effect of viscous dissipation [10] and compression work [11,12] on the thermal and hydrodynamic behavior is more significant and should not be neglected [13].
The kinetic-based LBM [14], due to its advantages, such as high computational efficiency, easy implementation, and use of parallel algorithms, has become a powerful numerical technique for simulating fluid flows and modeling the physics of fluids [15,16,17].Various thermal LB models or Boltzmann-based schemes have been employed to investigate natural convection problems [18,19,20].The multispeed LBM is a straightforward extension of the isothermal LBM; however, it usually suffers from severe numerical instability, and the simulation of the temperature variation is limited to a narrow range [21].In the hybrid LBM, the flow simulation is performed by using the LBM; however, the energy equation is solved through conventional numerical methods.Therefore, this method has not demonstrated the advantages of the standard LBM [19].The DDF-LBM uses two distribution functions, one for the velocity field and the other for the temperature or internal energy field.And it takes into account the viscous heat dissipation and compression work [22,23].However, most of the DDF-LBM models are decoupling models in the sense that the change of temperature cannot influence the velocity field [14].Hence, the decoupling models are only applicable to Boussinesq flows in which the temperature variation is minor.When these models are applied to the thermal problems in which the temperature field has significant effects on the flow field, the decoupling between the momentum and the energy transports causes considerable errors.
To overcome this problem, on the basis of the total energy DDF-LBM [24], we attempt to propose a coupled DDF thermal LBM (coupled DDF-TLBM) [25] that incorporates viscous heat dissipation and compression work.In order to realize the coupling between the momentum and the energy fields, this method introduces the temperature change into the LB momentum equation in the form of a momentum source, which affects the distribution of flow velocity and density.Although numerous studies have demonstrated the capabilities of LB models in simulating natural convection, an appropriate model for reflecting the effect of viscous heat dissipation and compression work is still desirable for microscale Rayleigh-Bénard convection without the Boussinesq approximation.

The Microscale Rayleigh-Bénard Convection
As shown in figure 1, the microscale Rayleigh-Bénard convection occurs in a 2D rectangular cavity with thermally insulated sidewalls filled with a viscous incompressible static fluid with an initial temperature T0.The temperatures of the top and bottom plates are constant but different.Specifically the top plate is cold with a temperature Tc, and the bottom plate is hot with a temperature Th.The height, width and aspect ratio of the rectangular cavity are denoted as H, L and A with A=H / L.
The Ra [26] and Prandlt number (Pr) [27] are dimensionless parameters that are associated with Rayleigh-Bénard convection flow.Ra is the ratio of buoyancy to viscosity forces multiplied by the ratio of momentum to thermal diffusivities.It characterizes the transition between the conduction-and convection-dominated flows.Pr is the ratio of viscous diffusion to thermal diffusion.
where cp is the specific heat coefficient at constant pressure, cp=(D+2)R/2 with D being the spatial dimension and R the gas constant, g is the acceleration of gravity, β is the isobaric coefficient of thermal expansion, µ is the dynamic viscosity, κ is the thermal conductivity, and ρ is the fluid density.The evolution equations can be written as follows: where x is the position of the lattice node; δt is the time step; 2 (2 ) , where f τ and h τ are the dimensionless relaxation times for momentum and total energy, respectively; i f and i h are respectively the density distribution function and total energy distribution function; eq i f and eq i h are the equilibrium distribution functions; i e and u are the discrete velocities of the fluid particles and macroscopic velocity, respectively; and i F and i q are two terms related to the external force: e u e a a u (5) 0 where a is the external force acceleration, 0, where ω is the weighting factor for the various lattice links, and 1 36 According to the conservation of the collision operator, the macroscopic variable velocity and total energy are directly calculated by using the distribution functions: , 2 2 Through a Chapman-Enskog analysis, it is found that the pressure p0 is the dynamic pressure rather than the thermodynamic pressure, and the change of the temperature field cannot influence the velocity field [24].
In order to realize the coupling between the momentum field and energy field, the coupled DDF-TLBM presented in this paper introduces the temperature change into the LB momentum equation in the form of the momentum source Ki, which affects the distribution of flow velocity and density.It is form of In accordance with the definition of the total energy 2 2 is the specific heat at constant volume, the temperature T at the lattice x and moment t is given by The temperature difference ( , ) i T t x between the lattice x and the surrounding lattice point in the direction i e is given by ( , ) ( , ) -( , ) The momentum source Ki produced by the temperature change is written as follows: e u e x x x (14) The temperature change in Ki includes two parts: one for the same lattice at different times that affects the thermal diffusion coefficient and viscosity coefficient of the fluid; and the other for different lattices at same time that affects the velocity distribution of the flow field.The temperature change is introduced into the LB equation in the form of the momentum source and influences the flow field through collision and migration.This coupled DDF-TLBM can realize the coupling between the momentum field and energy field, and overcome the limitation of the Boussinesq approximation.Hence it extends the application range of the thermal LBM.f τ and h τ are expressed by the dimensionless characteristic parameters Pr, Ra, and Mach number (Ma) ( In practical applications, the flow boundary conditions are usually specified in terms of the fluid variables.The nonequilibrium-extrapolation approach [24] is employed to transform thermos-hydrodynamic boundary conditions to the boundary conditions of the distribution functions due to its simplicity, second-order accuracy, and good robustness.In this method, the distribution function at a boundary node is separated into the equilibrium and nonequilibrium components.The equilibrium component is determined by using the macroscopic variables of the boundary, and the nonequilibrium component is identified by using the distribution function at the nearest node in the fluid region [14].According to this approach, the density and energy distribution function at a boundary node b  (18) where the internal equilibrium distribution functions 0 ( , )= ( , ) / 2

RESULTS AND DISCUSSION
By using the method discussed above, two models will be employed to describe the Rayleigh-Bénard convection.One considers the viscous heat dissipation and compression work (model I) and the other does not take into account the viscous heat dissipation and compression work (model II).The numerical-simulation experiments are conducted to investigate the effect of the viscous heat dissipation and compression work on the temperature distribution, flow distribution, and average Nu [28] for various values of the parameters, such as Ra and the aspect ratio A. In the computations 128×64 lattices is employed to divide the flow domain (128µm×64µm).Table 1 lists the default values of the parameters used in the simulations.As shown in figures 2 and 3, when Ra = 10 3 , the heat conduction plays a crucial role in the two models, and the convection effect is weak.The effect of the viscous heat dissipation and compression work is minor.Consequently, the velocity profiles and isotherms of the two models are similar.When Ra = 10 4 , the convection effect is enhanced, and the effect of the viscous heat dissipation and compression work becomes apparent.Therefore, the differences between the two models are evident.The velocity profile in figure 2 (a2) has larger vortexes and a clearer boundary than that in figure 2 (b2).Figure 3 (a2) and (b2) illustrates the differences in temperature distribution.When Ra = 10 5 , the differences between the two models are significant.In figure 2 (a3), the shape of the vortex is upside down.The convective heat transfer plays a crucial role in the model Ι.For the model Π, the role of heat conduction is the same as that of thermal convection.When Ra = 10 6 , the convection heat transfer plays a dominant role in both models; however, the model Ι has clearer velocity and temperature boundary layers.
This result shows that the effect of the viscous heat dissipation and compression work can promote convection heat transfer, which considerably affects the velocity and thermal distribution of microscale Rayleigh-Bénard convection.In heat transfer at a boundary (surface) within a fluid, Nu is the ratio of convective to conductive heat transfer across (normal to) the boundary.When the Nu is higher, the effect of the convective heat is more evident.To further explore the effect of viscous heat dissipation and compression work, Table 2 lists the average values of Nu ( Nu ) of the two models for various values of Ra, which are calculated as follows: As shown in Table 2, Nu increases with the increase of Ra.The total amount of heat transferring to the cold wall from the hot wall is increased; in other words, the convection resistance decreases, and the thermal conductivity resistance increases.For the same Ra number, Nu of the model Ι is higher than that of the model Π.This suggests that the viscous heat dissipation and compression work enhances the convective heat transfer.Therefore, the effect of viscous heat dissipation and compression work is crucial to Rayleigh-Bénard convection and should not be ignored.With the length fixed, the height of the flow domain is adjusted to obtain different aspect ratios (A). Figure 4 shows the velocity distribution profiles of the flow field at a steady state for A = 1, 1/2, 1/3, 1/4, 1/5, and 1/6. Figure 4 (a1)-(a6) is for the velocity distribution profiles of the mode Ι, and figure 4 (b1)-(b6) for the velocity distribution profiles of the model Π. Figure 5 displays the corresponding isotherm diagrams.
From figure 4 one can see that when A = 1, the velocity distributions and isotherms of the two models are similar.For the velocity field, a large eddy is formed in the square cavity.Note that the eddy of the model Ι is more intense than that of the model Π.The case for A = 1/2 has been discussed in Sect.3.1.The decrease in the system size enhances the effect of viscous heat dissipation and compression work.When A = 1/3, 1/4, 1/5, and 1/6, the differences between the two models are clearer.For the model Π, when A = 1/3, 1/4, 1/5, and 1/6 there are 3, 4, 5, and 6 eddies, respectively; in other words, the number of vortexes is inversely proportional to the aspect ratio A. By contrast, for the model Ι, when A = 1/3, 1/4, 1/5, and 1/6 there are 4, 6, 8, and 10 eddies, respectively.From figure 5, the similar phenomena are observed in the isothermal diagrams.As the aspect ratio decreases, the size of the flow domain decreases, and the effect of the viscous heat dissipation and compression work increases; the number of eddies increases linearly, and their shape is also changed.For the model Π, the eddies are elliptical, and the long axis of the eddies is tilted.For the model Ι, the long axis of the elliptical spiral is vertical.Therefore, the effect of viscous heat dissipation and compression work should be considered for the microscale flow.
The temperature profiles along the height of cavity for Ra=10 4 and A=1 are shown in figure 6.The blue lines are for the result obtained based on the model Ι and the red lines for the model Π. Figure 6 shows that the temperature gradient near the wall is large and the temperature change becomes smooth in the middle region, and the smooth region of model Ι is larger than that of model Π.Therefore, the convection effect is significant in the model Π due to the viscous heat dissipation and compression work.

CONCLUSIONS
In order to investigate the effect of viscous heat dissipation and compression work, the coupled DDF-TLBM is applied to two microscale Rayleigh-Bénard convection numerical models, one considering the viscous heat dissipation and compression work and the other not considering the viscous heat dissipation and compression work.The numericalsimulation experiments were conducted to study the effect of the viscous heat dissipation and compression work on the temperature distribution, flow distribution, and average Nu at different values of Ra and aspect ratios (A).It is found that when Ra > 104, the effect of viscous heat dissipation and compression work promote the convection heat transfer and have a profound effect on the Rayleigh-Bénard convection model; as the size of the flow domain decreases, the number of eddies increases linearly and their shape is also changed in the model that takes into account the viscous heat dissipation and compression work.
This study demonstrated that the coupled DDF-TLBM can provide a new sight for the microscale flow that incorporates viscous heat dissipation and compression work.Moreover, we demonstrated that viscous heat dissipation and compression work can promote the convection heat transfer and increase the number of vortexes, which is crucial for the microscale flow and therefore should be considered.These findings are crucial, and the model constructed in this study will be useful in promoting the development of microfluidic systems

Figure 1 :
Figure 1: Diagram of microscale Rayleigh-Bénard convection the microscale Rayleigh-Bénard convection model, T0 is the reference temperature ( 0 = ( ) 2 h c T T T + for the Rayleigh-Bénard convection model).In the following we use the D2Q9 lattice model.Thus the local fluid density and energy equilibrium distribution functions and the discrete velocities of the fluid particles i e are given by

2
constitute the coupled DDF-TLBM that incorporates viscous heat dissipation and compression work.For the flow model without the viscous heat dissipation and compression work, an internal energy distribution function i g is used in the temperature field, and the evolution equation is expressed as follows:

Figure 3
displays the isotherms corresponding to the velocity profiles.

Figure 6 :Figure 7
Figure 6: The normalized temperature profiles along the height of the cavity at Ra=10 4 and A=1, the blue line denotes the model Ι and the red line denotes the model Π

Figure 7 :
Figure 7: Relationship between the average Nu ( Nu ) and the aspect Ratio (A); the blue □ line denotes the model Ι, and the green ○ line denotes the model Π

Table 1 .
Default Values of the Basic Parameter Variables Used in the Simulations Impact of Ra on the Flow Field and Temperature Distribution

Table 2 :
Nu of the Two Numerical Models for Various Values of Ra Impact of Aspect Ratios (A) on the Flow Field and Temperature Distribution Temperature Distribution Next consider the effect of viscous heat dissipation and compression work on the flow field and temperature distribution of the microscale Rayleigh-Bénard convection with different geometries.