Numerical study of three-dimensional natural convection in a cubical cavity at high Rayleigh numbers

A systematic numerical study of three-dimensional natural convection of air in a differentially heated cubical cavity with Rayleigh number (Ra) up to 10 is performed by using the recently developed coupled discrete unified gas-kinetic scheme. It is found that temperature and velocity boundary layers are developed adjacent to the isothermal walls, and become thinner as Ra increases, while no apparent boundary layer appears near adiabatic walls. Also, the lateral adiabatic walls apparently suppress the convection in the cavity, however, the effect on overall heat transfer decreases with increasing Ra. Moreover, the detailed data of some specific important characteristic quantities is first presented for the cases of high Ra (up to 10). An exponential scaling law between the Nusselt number and Ra is also found for Ra from 10 to 10 for the first time, which is also consistent with the available numerical and experimental data at several specific values of Ra. 2017 Elsevier Ltd. All rights reserved.


Introduction
Natural convection flow (NCF) in a differentially heated cubical cavity is one of the fundamental flow configurations in heat transfer and fluid mechanics studies, and it has many significant applications, including air flow in buildings, cooling of electronic devices, and energy storage systems. In recent years, with the rapid advance of the computer technology, direct numerical simulation (DNS) has become a popular and competitive way to study thermal convection flow problems.
The early numerical studies of NCF were usually restricted to two-dimensional (2D) configuration with relatively low Rayleigh numbers ðRaÞ. The pioneering work of de Vahl Davis et al. [1] provided original benchmark solutions for a square 2D cavity with 10 3 6 Ra 6 10 6 ; afterward, more accurate results were presented by Hortmann et al. [2] using the multi-grid method with a much finer mesh. Many others have repeated results with Ra up to 10 8 [2][3][4][5][6].
As actual flow is always three-dimensional (3D), many efforts have also been made on 3D simulations. For example, Mallinson et al. [7] investigated the effects of a certain aspect of a ratio on flow patterns with Ra up to 10 6 ; Fusegi et al. [8] simulated the NCF in an air-filled cubical cavity for Ra of 10 4 and 10 6 , and clari-fied 3D structures of flow and temperature; Labrosse et al. [9] observed the hysteretic behavior by using a pseudo-spectral solver; the 3D cavity of aspect ratio 4 with periodic lateral walls was studied by Trias et al. [10,11], and significant differences were observed in flow dynamics between 2D and 3D results. They also emphasized that the NCF in a 3D cubical cavity with adiabatic lateral walls had received comparatively less attentions [8,12,9,13,14].
The above mentioned numerical simulations of NCF are performed by the traditional computational fluid dynamics (CFD) methods on the basis of the Navier-Stokes equations (NSEs), which are a set of second-order nonlinear partial differential equations (PDEs). Recently, kinetic methods based on the Boltzmann model equation have become an alternative method to the NSEs with some distinctive features. Different from the NSEs with a nonlinear and nonlocal convection term, the Boltzmann equation is a firstorder linear PDE, and the nonlinearity resides locally in its collision term. These features make kinetic methods easy to realize and parallelize with high computational efficiency. Many kinetic methods have been recently utilized to simulate NCF problem, such as the lattice Boltzmann methods (LBM) [15][16][17][18][19][20][21][22][23] and the gas-kinetic scheme [24].
However, up to date, the study of high Ra (up to 10 10 ) NCF in a differentially heated cubical cavity is limited to 2D configuration, while most of the available investigations on 3D NCF are for low ÀRa (up to 10 7  finer mesh adjacent to wall boundary to capture boundary layer than that in the center of cavity. Therefore, an accurate, stable and mesh flexible method is preferable for numerical study of the 3D NCF. Recently, starting from the Boltzmann model equation, a discrete unified gas-kinetic scheme (DUGKS) was proposed for both hydrodynamic and rarefied flows [25][26][27]. As a finite-volume (FV) method, DUGKS can be easily implemented on non-uniform or unstructured meshes to satisfy the local accuracy requirement [28][29][30]. Particularly, although sharing the common kinetic origin, some distinctive features also exist between DUGKS and LBM. In fact, several comparative studies of the standard LBM and DUGKS have been preformed systematically for laminar flows [28,31], turbulent flows [32,33], and natural convection flows [29] in previous work. Generally, for flows without solid boundaries, for example the decaying turbulent flow, the accuracy of standard LBM is slightly higher than the DUGKS [32], while for flows involving solid boundaries, the DUGKS is even more accurate than the standard LBM [28,33]. Furthermore, owing to the semi-implicitness in the construction of gas distribution function at the cell interfaces, the DUGKS is much more stable and robust than LBM [28,29,32]. However, with a same regular grid, the standard LBM is faster than the DUGKS per iteration [28,29]. But benefiting from the FV nature, non-uniform meshes can be easily employed without loss of accuracy and additional efforts in DUGKS, and its efficiency can be significantly improved by employing a non-uniform mesh according to the local accuracy requirement [28,29,32]. This is the main reason why we use the DUGKS, instead of the LBM, to study the high Rayleigh number natural convection flow, which requires much fine mesh near walls to resolve the thin boundary layers. Although the standard LBM has low numerical dissipation [34], it can only be implemented on regular meshes due to its special streaming process, and the existing reported studies of 3D NCF are limited to low Ra [17,23]. Some FV based LBM have been also developed in the past decades, but it has been demonstrated that the DUGKS is obviously superior to the current best FV-LBM [35] in terms of accuracy and numerical stability [31]. In addition, the kinetic nature makes DUGKS suitable for parallelized computing. High computational efficiency is essential for large scale 3D simulations. In order to simulate the incompressible thermal flow, a coupled DUGKS (CDUGKS) has been proposed using the double distribution strategy, and its accuracy, efficiency and numerical stability have been validated by simulating the 2D NCF with Ra up to 10 10 [29]. In this work, we will contribute to study the NCF in a differentially heated cubical cavity with Ra up to 10 10 using the CDUGKS. The method is firstly validated by comparing with available numerical and experimental data. Flow characteristics and heat transfer are then to be investigated. Finally, a scaling correlation between Rayleigh and Nusselt numbers with Ra up to 10 10 will be obtained for the first time.

Kinetic model equations
The coupled discrete unified gas-kinetic scheme is derived from the following Boltzmann model equations [24] @f @t @g @t þ n Á r x g ¼ W where f and g are gas distribution functions for velocity and temperature fields, respectively, and both are functions of space x, time t, and molecular velocity n; f eq and g eq are the corresponding equilibrium states here R is the gas constant, T 1 and T 2 are the constant variances. For convenience, we set T 1 ¼ T 2 in this study. s v ¼ m=RT 1 and s c ¼ j=RT 2 are the corresponding relaxation times, here m and j are, respectively, the kinematic viscosity and heat conduction coefficient, which determine the Prandtl number Pr ¼ m=j. For low speed flows, the external force term F can be approximated as [36] F ¼ Àa Á r n f % Àa Á r n f eq ¼ here a is the acceleration due to buoyancy force, and is approximated by the Boussinesq assumption where g 0 is the gravitational constant,ĝ is the unit vector in the gravitational direction.
In computation, the continuous molecular velocity space should be approximated by a discrete velocity set fn i ji 2 Zg, so that the integration on molecular velocity space can be numerically computed. For nearly incompressible flow (i.e., when the Mach number Ma ( 1), the equilibrium states can be approximated using the Taylor expansion to the second order, i.e., where f eq i ¼ x i f eq ðn i Þ; g eq i ¼ x i g eq ðn i Þ; and W i is the weight coefficient corresponding to molecular velocity n i . In the present study, we use nineteen velocities in three dimensions, i.e., the D3Q19 model, with where c ¼ ffiffiffiffiffiffiffiffiffiffiffi 3RT 1 p , and the corresponding weight coefficients are W 0 ¼ 1=3; W 1;...;6 ¼ 1=18 and W 7;...;18 ¼ 1=36. The discrete distribution functions f i ðx; tÞ ¼ x i f ðx; n i ; tÞ and g i ðx; tÞ ¼ x i gðx; n i ; tÞ satisfy the following equations The fluid density, velocity, and temperature can be obtained from the discrete distribution functions,

DUGKS for velocity field
The DUGKS is a FV method in which the computational domain is divided into a set of control volumes. We integrate Eq. (10) over a control volume V j centered at x j from t n to t nþ1 (the time step Dt ¼ t nþ1 À t n is assumed to be a constant), and use the midpoint and trapezoidal rules for time integrations of the convection and collision terms, respectively, so we get the updating equation as [25] f nþ1 is the micro-flux across the cell interface, and Since the collision term X i is mass and momentum conservative, the density q and velocity u can be computed by The key ingredient in updatingf i is to evaluate the interface flux F nþ1=2 i , which is solely determined by the distribution function f i ðx; t nþ1=2 Þ there. In DUGKS, after integrating Eq. (10) along the characteristic line of molecule within a half time step ðh ¼ Dt=2Þ, the evaluation of f i ðx; t nþ1=2 Þ can be traced back to the interior of neighboring cells, where ; t n Þcan be approximated by linear interpolation.
Again using the conservative property of X i , the density q and velocity u at the cell interface can be obtained from from which the equilibrium distribution function f eq i x b ; t n þ h ð Þis determined. Therefore, based on Eq. (18) and the obtained equilibrium state, the original distribution function f i at the cell interface can be determined from from which the micro-flux Eq. (14) can be evaluated.
In computation, we only need to follow the evolution off i in Eq. (13). The required variables for its evolution are determined by 2.3. DUGKS for temperature field DUGKS for Eq. (11) can be constructed similarly as for Eq. (10). Eq. (11) is firstly integrated in the same control volume V j from t n to t nþ1 , then, the same integration rules are employed to approxi-mate the time integrations of convection and collision terms, one can get the updating equation for g i , where F is the micro-flux, are auxiliary distribution functions. Based on the conservative property of W i , the temperature can be computed as In order to evaluate the micro-flux F, we again integrate Eq. (11) within a half time step h along the characteristic line with an end point x b at the cell interface, and use the trapezoidal rule to approximate the time integration of collision term, and Þis made around the cell interface x b . Based on Eq. (28) and the compatibility condition, the temperature at the cell interface can be computed as Together with the conserved variables in velocity field, the equilibrium distribution function g eq i ðx b ; t n þ hÞ can be fully determined. Then, the original distribution function g i can be obtained, from which the interface numerical flux F is determined. In computation, it only needs to follow the evolution ofg according to Eq. (23). The required variables for its evolution are determined by

Kinetic boundary conditions
In this study, kinetic boundary conditions for velocity and temperature fields are required. For the velocity field, no-slip boundary on the fixed wall (x w ) can be realized by the bounce-back rule, which gives where n is the unit vector normal to the wall pointing to the cell. As for the temperature field, two types of boundaries i.e., fixed temperature and adiabatic walls are considered. For the constant temperature boundary, the distribution function for molecule leaving the wall can be given by [37] where T w is the wall temperature. The adiabatic boundary, which is a Neumann boundary condition, can be realized by the bounceback rule [21], The distribution function g i ðx w ; t þ hÞ for molecules moving towards the wall, i.e., n i Á n 0, can be obtained following the procedure described in Section 2.3.

Algorithm
For clarify, we now list the computational procedure of DUGKS.
Specifically, with initializedf 0 j;i andg 0 j;i in all cells centered at x j , the procedure of the DUGKS at each time t n reads as follows:

Numerical results
In this work, natural convection of air in a differentially heated cubical cavity with the adiabatic lateral walls is studied in 3D. As illustrated in Fig. 1, the configuration is a cubical box with a cold wall ðT ¼ T c Þ on the left side (x ¼ 0) and a hot wall ðT ¼ T h Þ on the right side (x ¼ H), and the other four walls are adiabatic. The gravity is along the yÀdirection. The Rayleigh number is defined as Ra ¼ g 0 bDTH 3 =ðmjÞ, where DT ¼ T h À T c is the temperature difference between the hot and cold walls, and H is the length of the cavity.
In the simulations, we set the length of the box H ¼ 1, the Prandtl number Pr ¼ 0:71, the temperature T c ¼ 0 and T h ¼ 1 on the cold and hot walls, respectively. In addition, we set RT 1 ¼ 10 is small enough to satisfy the incompressible approximation, where is the characteristic velocity. The velocities are normalized by the characteristic velocity u 0 unless otherwise stated. The velocity boundary condition, Eq. (33), is applied to all walls; the temperature boundary conditions, Eqs. (34) and (35), are applied to the isothermal and adiabatic walls, respectively. The time step is determined by the Courant-Friedrichs-Lewy (CFL) condition, i.e., Dt ¼ gDx min =n max , where g is the CFL number, Dx min is minimum grid spacing, and n max is the maximum discrete molecular velocity. We set the CFL number g ¼ 0:95 in the following simulations.
The FV nature of DUGKS make it easy to vary mesh resolution according to the local accuracy requirement. In our simulations, a set of non-uniform meshes are adopted with N grid points in each direction, and the mesh gradually becomes finer from the center to walls of the cavity. The location of a control volume center in which a is a constant that determines the mesh distribution. Here a is set to 2:5, and the mesh distribution is shown in Fig. 2. In our simulations, a suitable mesh should be adopted to make a balance between computational accuracy and efficiency. For the cases of In order to make a quantitative study, we will measure some important quantities of the NCF, including the maximum horizontal velocity at the vertical centerline of the midplane (z ¼ H=2), u max , and the corresponding yÀcoordinate, the maximum vertical velocity at the horizontal centerline of the midplane (z ¼ H=2), v max , and the corresponding xÀcoordinate, the mean Nusselt num-  Nu l ðy; zÞj x¼0 or x¼H dy; Nu l ðy; zÞ ¼ @Tðy; zÞ @x ; ð37Þ in which Nu l is the local Nusselt number.

Validation
The comparative study with the available results in the literature is firstly performed to validate the 3D CDUGKS code. Initially, we set the flow velocity u ¼ 0 and the temperature T ¼ 0:5 for the whole flow field. To assess whether the flows have reached the steady-state, the following criterion is applied, ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P uðtÞ À uðt À 1000DtÞ k k 2 q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi P kuðtÞk 2 q < 10 À8 : ð39Þ Table 1 shows the quantities of interest on the symmetry plane of z ¼ 0:5, as well as the mean Nu (Nu m ) and the overall Nu (Nu o ) on the hot wall for 10 3 6 Ra 6 10 6 . The results of a high-resolution finite-difference (FD) method [8] are included for comparison. As can be seen, good agreements are achieved with the benchmark results, and the maximum relative error is less than 5%, which is acceptable for engineering applications. Also, it is found that similar to the LBM [17], the CDUGKS almost overall underestimates the results of these quantities when compared with the same reference solution. The deviation can be attributed to the difference in the numerical methods. The reference solutions were obtained by a third-order QUICK scheme for the incompressible Navier-Stokes equations on a fine non-uniform mesh of 62 3 , while the DUGKS is a second-order accurate method, and a relatively coarse nonuniform mesh of 50 3 is employed in the validation. However, the deviations between the results of DUGKS and reference solutions are still acceptable (less than 5%).
In addition, we compare the temperature and velocity profiles given by CDUGKS with those from the experiments or the NS solver that are available at some specific values of Ra. It should be noted that cavities with large depth aspect ratios were usually employed in the experiments, which precludes the precise comparisons of each set of data. Fig. 3 shows the temperature profiles on the symmetry plane z ¼ 0:5 for the case of Ra ¼ 10 5 , in which the experimental data at Ra ¼ 1:89 Â 10 5 acquired by a Mach-Zehnder interferometry technique [38] and the NS solver results obtained by FD method [8] are also included for comparison. It is found that Fig. 5. Isotherms on the symmetry planes for 10 3 6 Ra 6 10 6 . the respective temperature distributions at the mid-lines y ¼ 0:5 and x ¼ 0:5 agree well with the FD results [8], and are in good agreement with the experimental data at the mid-line y ¼ 0:5. However, there is a noticeable discrepancy at the mid-line x ¼ 0:5 between the two numerical and experimental results, and the derivation increases when approaching the horizontal walls. This may be attributed to the unavoidable heat transfer through adiabatic walls in the experimental situations [38] and slight difference in Ra between the numerical simulations and experiment. Fig. 4 shows the velocity profiles on the symmetry plane z ¼ 0:5 for the case of Ra ¼ 10 5 . In addition to the above mentioned reference results, another set of experimental data at Ra ¼ 1:03 Â 10 5 [39] are also included. It is clearly seen that the CDUGKS results are in excellent agreement with the FD results, and agree reasonably well with two sets of experimental data, despite slight deviations in the peak values of velocities.
3.2. Flow characteristics and heat transfer for 10 3 6 Ra 6 10 10 The flow characteristics and heat transfer are two important features of the 3D NCF, but systematic investigation is still lack in the literature, especially for the cases of high Ra (larger than 10 6 ). In this subsection, a detailed study of the flow characteristics and heat transfer of NCF in a 3D cavity is performed for 10 3 6 Ra 6 10 10 .
Figs. 5 and 6 show isotherms or instantaneous isotherms on the symmetry planes for 10 3 6 Ra 6 10 6 and 10 7 6 Ra 6 10 10 , respectively. It can be seen that as Ra increases, isotherms change from almost vertical to be horizontal in the center of the symmetry plane z ¼ 0:5, and are vertical only adjacent to the isothermal walls, which indicates that the dominant heat transfer mechanism changes from conduction to convection. In addition, temperature boundary layers are developed near the isothermal walls, and become thinner as Ra increases. However, such phenomenon does not appear near the adiabatic walls. Furthermore, as Ra reaches to 10 9 and 10 10 , the isotherms become oscillatory, which means that the heat transfer has turned to be time-dependent. Moreover, as shown on the symmetry planes of y ¼ 0:5 and x ¼ 0:5, the temperature adjacent to the lateral adiabatic walls is lower than that in the center region, suggesting that the adiabatic walls inhibit heat transfer, but such effects reduce as Ra increases. It is also noted that for the cases of Ra 6 10 6 , isotherms on the symmetry plane of z ¼ 0:5 are very similar to the 2D results [17,40]. nels are developed near the center of walls, and as Ra increases, they move to the corners of the cavity. For example, as shown in Fig. 8ðaÞ, the maximum velocity is located adjacent to two corners of the isothermal walls, which means two tunnels have been developed there. In addition, we also find that similar to the temperature field, velocity boundary layers are formed near the isothermal walls, which become thinner with increasing Ra, but no apparent velocity boundary layers are developed adjacent to the adiabatic walls. Furthermore, when the Ra reaches 10 9 , the flow field is evolving to chaos. Heat transfer in the cubical cavity is characterized by the Nusselt number (Nu), which is a measure of the relative strengths of the convective and conductive heat transfer processes. The heat transfer is dominated by conduction in the flow domain where Nu < 1, and by convection as Nu > 1. Fig. 9 shows the contours of the local Nusselt number (Nu l ) on the hot wall. It can be seen that the value of Nu at the symmetry plane of z ¼ 0:5 is higher than the other regions, suggesting that the lateral adiabatic walls suppress the convective heat transfer in the 3D cavity, but such effect is gradually reduced as Ra increases. It is also observed that the higher local Nusselt number appears near the bottom of hot wall, which means that convection is the dominant heat transfer process in this region. This is because that with anticlockwise vortex rolling, the cold fluid near the bottom of cavity is closer to the hot wall than that in the vicinity of upper region, see Figs. 5 and 6, so the convective heat transfer is more intensive at the bottom of hot wall.
Quantitatively, we present time-averaged characteristic quantities of interest on the symmetry plane of z ¼ 0:5 for 10 7 6 Ra 6 10 10 in Table 2, which is absent from the literatures. Note that the provided results are already mesh independent, and the time-averaged operation is preformed over a long enough period. It is found that the maximum horizontal velocity is located near the bottom adiabatic wall, and the maximum vertical velocity moves to the hot wall as Ra increases. Also, the vertical velocity near the isothermal hot wall is about one order of magnitude larger than the horizontal one near the adiabatic wall, which indicates that contrary to the results near the adiabatic wall, a boundary layer is formed adjacent to the isothermal wall. This result is consistent with the phenomenon we have observed in Fig. 7. Moreover, the conservation of heat transfer in the cubical cavity is critical to the DNS of high Ra NCF. Therefore, the mean and overall Nusselt numbers on cold and hot walls are included in Table 2 to validate the conservation property of our simulation. It is found that the maximum difference of Nu m or Nu 0 between the cold and hot walls is less than 0:5%, suggesting that the heat transfer conservation in the cavity is well predicted by the CDUGKS.

Nu-Ra correlation
Ra and Nu are two important parameters in description of heat transfer. For the 3D NCF in a cubical cavity, it has been demonstrated that there is an exponential relation between Nu and Ra In Fig. 10, the fitting relationships of Nu m À Ra (Eq. (40)) and Nu o À Ra (Eq. (41)), as well as the numerical results computed by CDUGKS, are shown. The available numerical and experimental data at some specific values of Ra are also included. As shown, the predicted results given by Eqs. (40) and (41) are in good agreement with the solutions obtained by the spectral Chebyshev method [13,14] and the experiment data [14], which validates the proposed scaling laws. Note that the scaling relationship for 10 7 6 Ra 6 10 10 reported here is given for the first time, which may help in thermal engineering applications.

Conclusions
In this paper, the natural convection of air in a differentially heated cubical cavity with the adiabatic lateral walls is studied by the coupled discrete unified gas-kinetic scheme (CDUGKS). The validation of our 3D code is carefully performed. It is shown that the CDUGKS results are in good agreement with those from the traditional CFD solvers and experiments. Particularly, the CDUGKS can accurately capture the temperature and velocity boundary layers adjacent to the isothermal walls. In addition, some important quantities computed by the CDUGKS also agree quantitatively well with the available numerical and experimental data.
The flow characteristics are studied systematically for 10 3 6 Ra 6 10 10 . It is found that the temperature and velocity boundary layers are formed near the isothermal walls, and become thinner as Ra increases, while no apparent boundary layer is developed near the adiabatic walls. In addition, it is observed that flow tunnels are developed close to the center of isothermal walls, and they move to the cavity corners in the joint of isothermal and adiabatic walls as Ra increases, meanwhile, the flow evolves from the steady-state to time-dependent state. When Ra approaches to 10 9 , the flow motions turn to be fully turbulent. Moreover, it is also  found that the lateral adiabatic walls have an inhibition on the distribution of temperature in the cavity. By studying the local Nusselt number distribution on the hot wall, we observe that the lateral adiabatic walls apparently suppress the heat transfer. Particularly, the Nusselt number on the hot wall has an instant drop near the adiabatic walls, which means that the convective intensity in the center of cavity is stronger than that close to the adiabatic walls. But, the effect of adiabatic walls on overall heat transfer decreases as Ra increases. Moreover, the time-averaged characteristic quantities of interest on the symmetry plane for the Rayleigh number up to 10 10 is presented for the first time. The results of Nusselt number on the cold and hot walls show that the heat transfer predicted by our simulations is well conserved in the cavity.
Finally, the scaling laws between the Nusselt number (local and overall) and Rayleigh number in the air filled cubical cavity for Ra up to 10 10 are revealed. The predicted results agree well with the available numerical and experimental data at several specific  Table 2 Time-averaged characteristic quantities on the symmetry plane of z ¼ 0:5 for 10 7 6 Ra 6 10 10 . Here Numcand Nu mh denote the mean Nu on the cold and hot walls, respectively. Likewise, Nuoc and Nu oh denote the overall Nu on the cold and hot walls, respectively. The velocities are normalized by the characteristic velocity u0 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gbDTH p .