Abstract
Different numerical methods have been implemented to simulate internal natural convection heat transfer and also to identify the most accurate and efficient one. A laterally heated square enclosure, filled with air, was studied. A FORTRAN code based on the lattice Boltzmann method (LBM) was developed for this purpose. The finite difference method was applied to discretize the LBM equations. Furthermore, for comparison purpose, the commercially available CFD package FLUENT, which uses finite volume Method (FVM), was also used to simulate the same problem. Different discretization schemes, being the first order upwind, second order upwind, power law, and QUICK, were used with the finite volume solver where the SIMPLE and SIMPLEC algorithms linked the velocity-pressure terms. The results were also compared with existing experimental and numerical data. It was observed that the finite volume method requires less CPU usage time and yields more accurate results compared to the LBM. It has been noted that the 1st order upwind/SIMPLEC combination converges comparatively quickly with a very high accuracy especially at the boundaries. Interestingly, all variants of FVM discretization/pressure-velocity linking methods lead to almost the same number of iterations to converge but higher-order schemes ask for longer iterations.
1. Introduction
Studying heat transfer and fluid flow using computational methods is easier [1], safer [2], and much less costly [3] compared to experimental techniques. There are a large number of problems which can be simulated with great accuracy to replicate experiments with high resolutions [4]. There are currently a range of approaches with the potential to serve in modeling heat transfer and fluid flows, such as the finite difference method (FDM), finite element method (FEM), finite volume method (FVM), lattice boltzmann method (LBM), boundary elements method (BEM), molecular dynamics simulation, and direct simulation Monte Carlo. The most widely employed approaches in the field of thermofluids are the first four [5]. However, application of FDM can be difficult when complex geometries are involved [6]. The FEM schemes can be intricate for solving conservative equations, while the nonstandard FEMs have low computational efficiency [7]. Application of FVM is difficult and complex to cases with complex moving boundaries [8]. LBM is a compressible model for ideal gases and can theoretically always simulate the compressible Navier-Stokes equations. With the Chapman-Enskog expansion [9], LBM can simulate incompressible flow for low Mach numbers (Ma < 0.15) albeit at the expense of a compressibility error [10, 11]. Besides, regular square grids used with LBM make it very hard to extend the simulation to curved boundaries [12]. All in all, the accuracy of all these numerical approaches is dependent on the problem configuration, discretization scheme, and numerical algorithm used [5]. As such, an important question to answer is about finding the best approach to solve a certain problem subject to computational efficiency and accuracy as the most important constrains. Along these lines, Rouboa and Monteiro [13] investigated the heat transfer phenomenon during cast solidification in a complicated configuration by FVM and FDM. A comparison between the numerical results and experimental ones indicated that both discretization approaches produced good outcome, with FVM being slightly better as it uses more information than FDM to capture spatial temperature variations. Despite recent progress in computing power and techniques, the literature review indicates a lack of comprehensive studies on selecting the ideal means of analyzing internal heat transfer and fluid flow problems. In particular, an optimal solution technique and procedure to simulate internal natural convection are yet to be presented. To fill this gap in the literature, laminar natural convection heat transfer of air inside a laterally heated square enclosure is investigated using both FVM and LBM. The simulation results were compared against those from the literature. Particular attention was given to different discretization techniques as well as pressure-velocity linking approaches to find the best method for simulating internal free convection problems.
2. Governing Equations
2.1. Finite Volume Method
Continuity, momentum, and energy equations were employed for flow analysis in a system depicted by Figure 1. Density was computed by invoking the Boussinesq approximation for [14]. The governing equations are written as follows [15].
Continuity equation: Momentum equations in and directions: Energy equation:
2.2. Lattice Boltzmann Method
The hydrodynamic and thermal Boltzmann equations with using density-momentum and internal energy distribution functions (double population) are as follows [16, 17]: Double population LBM model (TLBM) uses two separated distribution functions and for hydrodynamic and thermal fields, respectively. This model is the latest one among different presented models of thermal LBMs. In addition, it shows more accuracy and stability during the solution process. As LBM solution process naturally tends to divergence having a stable approach like TLBM helps the convergence. Microscopic velocities for a D2Q9 lattice model are [12] Heat dissipation and hydrodynamic and thermal equilibrium distribution functions are given by where and is the weight function. Equation (4) in discretized forms [18] reads The two last equations are implicit. Thus, the new functions and are developed to address this problem: Collision and streaming steps of LBM are simulated by applying (7)–(9) as follows: Finally, the hydrodynamic and thermal variables can be obtained as
3. Boundary Conditions
Figure 1 illustrates a schematic of the configuration analyzed in the present study along with the boundary conditions.
The nonequilibrium bounce-back model is used to simulate the no-slip boundary condition on the walls in LBM. This model improves accuracy compared to the usual bounce-back boundary condition and satisfies the zero mass flow rates at nodes on the wall. The collision occurs on the nodes located at the solid-fluid boundaries and distribution functions are reflected in a suitable direction, satisfying the equilibrium conditions [19].
The macroscopic boundary conditions for the present study are
4. Numerical Procedure
Our FVM solver uses the implicit line-by-line tridiagonal matrix algorithm [20, 21] to linearize the system of algebraic equations. First order upwind [22], second order upwind [23], power law [24], and Quadratic Upstream Interpolation for Convective Kinetics (QUICK) [25] schemes were applied in different trails to solve the same problem while the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) [26, 27] and SIMPLE-Consistent (SIMPLEC) [28, 29] procedures were selected for pressure-velocity coupling. The convergence criterion, maximum absolute error in each dependent variable, was set at 10−7.
In LBM, the zero values for , , and are applied as the initial conditions. However, to avoid problems in estimating the macroscopic variables in (12), the initial fluid density is set to unity. LBM dimensionless numbers Re, Ra, and Pr are defined identical to those of classical Navier-Stokes equations. However, the macroscopic numerical value should be calculated beforehand. For example, for Pr, one has the kinematics viscosity and thermal diffusivity determined in LBM as and , where and are hydrodynamic and thermal relaxation times and is the gas constant. The Prandtl number can then be written as . For the values of and are now known based on relaxation times, while the numerical values of , , , are predetermined and fixed.
4.1. Gravity Effects in LBM
The Boussinesq approximation was used as to give buoyancy force per unit mass defined as and . Hence, the discretized Boltzmann equation is written as Applying (8) and taking into consideration the effects of gravity, one has while for thermal macroscopic variables (11) is applied.
4.2. Deriving Navier-Stokes Equations from LBM
In order to derive Navier-Stokes equations from the incompressible lattice Boltzmann equation by using Chapman-Enskog expansion the discretized form of Boltzmann equation can be written as With as a small (perturbation) variable, the Chapman-Enskog expansion for and reads None of the nonequilibrium parts of the above equations should be used for estimating the macroscopic properties and : Using these equations together with Tailor expansion of Boltzmann equation around , the terms which are smaller than () dropped, and then substitute into (15), we have Macroscopic density and velocity variables can be achieved by applying the first and second order of moments, leading to where Amount of is determined by using and then using the zero and first order of moments of (18) together with : where Finally, making use of , at incompressible limit and ignoring the term in (21), continuity and momentum equations are recovered. In addition, the thermal energy equation would be recovered in a similar way; see [30, 31] for more details.
5. Grid Independence
Structured nonuniform grid distributions were applied for FVM simulations with a grid cluster near the walls to capture sharp velocity and temperature gradients. For LBM simulations structured grids based on D2Q9 lattice are applied. Extensive grid independence checks were performed, as indicated by Tables 1 and 2, to observe that a grid with 47961 and 67600 cells for all FVM solvers and LBM, respectively, leads to mesh-independent results.
6. Results
Our numerical results, from different solvers, were compared with benchmark experimental data from Krane and Jessee [32] as well as the numerical predictions of Khanafer et al. [34], Oztop and Abu-Nada [35], and Bakhshan and Emrani [33]. The main dimensionless parameters were the Rayleigh and Prandtl numbers, which are constant at and 0.71, respectively. The fluid thermophysical properties, as well as dimensionless numbers, are shown in Table 3.
For and , the dimensionless temperature and vertical velocity profiles at midheight are plotted in Figures 2 and 3 and contrasted with the results from [32–35]. These figures illustrate a superior adaptation between the present simulation results using the FVM and LBM models and those of [32–35] works. Although previous research shows that, for complicated turbulent fluid flow problems, the QUICK/SIMPLEC is the most accurate choice [33], Figures 2(a) and 3(a) indicate that for laminar internal convection heat transfer problems there is no dramatic difference among the studied discretization approaches. However, it is obvious from Figures 2 and 3 that the FVM results are more accurate than those of LBM. This could be attributed to the compressible nature of LBM [36, 37], which creates a compressibility error for incompressible flows [12]. Among the discretization/pressure-velocity linking approaches examined, 1st order upwind/SIMPLEC has the closest results to experimental benchmark data, especially for the temperature contours in the range . With vertical velocity distribution, however, the difference among FVM approaches is quite negligible. Nevertheless, the fact that the accuracy and stability of the convective terms comprise a contrasting pair is a general perception in the field of computational heat transfer. For instance, the first order upwind scheme is entirely stable even with strong false diffusion [38], while the second or third order schemes like QUICK are conditionally stable [25].
| (a) |
| (b) |
| (c) |
| (d) |
| (e) |
| (a) |
| (b) |
| (c) |
| (d) |
| (e) |
Table 4 successfully compares our numerical results with those available in the literature under similar conditions and geometry over a range of Ra values with . Slight discrepancies are observed in this table between some of the present work results and those of [34, 39–42] because of the differences between the employed discretization methods, as well as mesh generation types, as one would expect.
Table 5 provides the comparison of number of iterations and required CPU usage time for the different discretization methods considered here. As seen, LBM may take 4-5 times longer to converge and 8-9 times more iterations compared to FVM. There are two reasons for this. The first one is attributed to the way LBM handles heat transfer. Although in the present work the appropriate internal energy distribution function, , [43] was used to obtain the temperature field, this model even tends to diverge. Furthermore, with LBM modeling the corners ask for a large number of fine grids near the corners. These two matters cause the LBM solutions to be comparatively more time consuming.
According to Table 5, the number of iterations for all FVM discretization method/pressure-velocity linking approaches is nearly equal. In this case, the difference between the QUICK/SIMPLEC method that necessitates the largest number of iterations and the lowest one (power law/SIMPLE) is only 79 iterations, that is, a 5.2% difference. With respect to CPU usage time, these proportions are to some extent different. For example, when comparing the most time consuming method (QUICK/SIMPLEC) with the 1st order upwind/SIMPLEC approach, this time disparity is about 4.94%, while the number of iterations differs by only 1.65%. As expected, higher-order accurate schemes are more time consuming.
The effects of the solution method, discretization scheme, and pressure/velocity coupling approach on the streamlines and -velocity are illustrated by Figures 4 and 5. Two elliptical vortexes generally appear at the center of the cavity as a predominant feature of buoyancy-induced flow in a laterally heated square enclosure. In this context, the 1st order upwind scheme has the most precise results among the studied discretization schemes, especially in the northwest and southeast sides of the enclosure. Regarding the pressure/velocity coupling approaches, the maximum stream function values for the SIMPLE and SIMPLEC approaches are 0.0003 and 0.0003066, respectively, translating into 2.2% difference while the CPU usage time difference is only 4 s. For LBM, the value of stream function is 0.000278.
(a) 1st order upwind/SIMPLE
(b) 2nd order upwind/SIMPLE
(c) 1st order upwind/SIMPLEC
(d) 2nd order upwind/SIMPLEC
(e) Power law/SIMPLE
(f) QUICK/SIMPLE
(g) Power law/SIMPLEC
(h) QUICK/SIMPLEC
(i) LBM
(a) 1st order upwind/SIMPLE
(b) 2nd order upwind/SIMPLE
(c) 1st order upwind/SIMPLEC
(d) 2nd order upwind/SIMPLEC
(e) Power law/SIMPLE
(f) QUICK/SIMPLE
(g) Power law/SIMPLEC
(h) QUICK/SIMPLEC
(i) LBM
For the velocity contours in the direction Figure 5 shows that the -velocity contours have cross-diagonal similarity towards the axis. Thus, all the methods analyzed present comparable results with no obvious difference.
Figure 6 demonstrates the local Nusselt number distributions along the left hot wall. For all discretization schemes, the Nusselt number is high near the bottom of the left wall (because of extreme temperature variations) and declines towards the top of the wall. The comparison between different solvers reveals that the 1st order upwind scheme predicts the maximum Nusselt number while LBM leads to the lowest one with some fluctuations along the hot wall. Interestingly, LBM uses about 40% more grids in that region compared to FVM ones.
| (a) |
| (b) |
| (c) |
| (d) |
| (e) |
7. Conclusions
Numerical tests using the finite volume and lattice Boltzmann methods with various discretization schemes and pressure-velocity linking algorithms were conducted to obtain the optimum discretization/linking approaches to address the internal convective heat transfer problems. The flow and temperature fields, as well as number of iterations and solving time, were evaluated.
The significant observations made in this study are summarized as follows.(1)The finite volume method results are more accurate compared to those of LBM, especially at the corners.(2)LBM needs a 4-5-fold CPU usage time and 8-9 times more iterations compared to the finite volume method to solve the problem considered here.(3)Among the studied discretization/pressure-velocity linking algorithms, the 1st order upwind/SIMPLEC provides the most precise results against experimental benchmark data, especially in the boundary layers.(4)The numbers of iterations for all FVM discretization/pressure-velocity linking methods are nearly equal.(5)The higher-order accurate schemes are more time consuming.One, however, notes that the above observations are valid within the limits of the parameters and problem considered in this study and could not be generalized to other cases without further investigations.
Nomenclature
| , | Cartesian coordinates (m) |
| Density-momentum distribution function | |
| , | Enclosure height and width (m) |
| Gr | Grashof number |
| Gravitational acceleration (m s2) | |
| Heat dissipation | |
| Internal energy distribution function | |
| Microscopic velocity vector | |
| Pr | Prandtl number ( ) |
| Pressure (N m2) | |
| Ra | Rayleigh number (Gr Pr) |
| Specific heat capacity (J kg1 K1) | |
| Temperature, time (K), (S) | |
| Thermal conductivity | |
| Velocities vector and its components in and directions (m s1). |
| Dynamic viscosity (Pa S) | |
| Density (kg m3) | |
| Internal energy relaxation times | |
| Kinematics viscosity (m2 s1) | |
| Momentum relaxation times | |
| Thermal expansion coefficient (K1) | |
| Thermal diffusivity, - direction components (m2 s1). |
| Cold wall | |
| Equilibrium distribution function | |
| Hot wall | |
| Lattice velocity direction. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors gratefully acknowledge the High Impact Research Grant UM.C/HIR/MOHE/ENG/23, UMRG Grant RP012C-13AET and the University of Malaya, Malaysia, for support in conducting this research work.