Using the Lattice Boltzmann Method for the numerical study of non-fourier conduction with variable thermal conductivity

1. Introduction

 In many different fields of applied sciences and engineering, heat sources such as lasers and microwaves have been extensively used under extremely short durations or very high frequencies. In such situations, the classical Fourier’s heat conduction theory becomes an inaccurate approach to describing diffusion and predicting temperature distribution. The non-Fourier effect is a more reliable means through which to explain the aforementioned [.1] processes. Non-Fourier motion (visually perceived motion that cannot be explained simply on the basis of the auto correlational structure of a visual stimulus) is a well-recognized phenomenon that is generally considered to be due to the non-linear preprocessing of a visual stimulus prior to standard motion analysis. In 1822, Fourier first investigeted an experimental relation between conduction heat transfer rate of a substance and the temperature gradient in the direction of energy flow. He deduced that the heat flux resulting from conduction heat transfer is proportional to the magnitude of the temperature gradient and opposite to it in sign[.2]  [1]. This observation may be expressed for unidirectional conduction heat transfer as Eq. (1):

        Most thermal conduction or mass diffusion problems are depicted and investigated using Fourier’s law, on whose basis non-physical conclusions are derived regarding situations that involve very high temperature gradients, extremely short durations, very low temperatures, and very small structural dimensions. The use of Fourier’s law is complemented with a modified flux model of transfer processes with finite speed waves[.3] . For an improved explanation of heat conduction in solids, however, researchers were driven to develop non-Fourier heat conduction theories. The mathematical statement of non-Fourier heat conduction law, which represents the time lag of heat waves, is a type of differential equation.

        If the occurrence of temperature difference  across depths  and  is the cause of heat conduction, the propagation of energy via conduction at any locations x and y in a medium with thermal conductivity k is expressed thus[.4] :

where  is the thermal lag time, that is, the rate of thermal diffusivity ( ) to the square of the speed of heat wave propagation , .

        The non-Fourier heat conduction equation has found broad applications in the analysis and design of thermal systems [1–3]. recently, the lattice Boltzmann method has become a powerful numerical algorithm for simulating heat transfer phenomena and fluid flows. The LBM is based on simplified mesoscopic equation engineering and can easily incorporate underlying physics into numerical solutions. The method is a powerful technique for the computational modeling of various complex fluid flow and heat transfer problems in complex geometries. It is also a discrete computational approach grounded in the Boltzmann equation. It regards a typical volume element of fluid as composed of a collection of particles that are represented by a particle velocity distribution function at each grid point. Time is counted in discrete time steps, and particles are assumed to collide with one another as they move, possibly under applied forces.

        Despite the advantages of the non-Fourier heat conduction mechanism, there are few studies about it. Ho et al. [4] and Mishra et al. [5, 6], investigated non-Fourier heat conduction in one-dimensional planar geometry. In some studies also analyzed non-Fourier heat conduction in planar and rectangular geometries [7–11] and cylindrical and spherical geometries [12–15]. Variable thermal conductivity heat transfer occurs in many engineering applications. Such variability can be observed, for example, during heat transfer in furnaces, boilers, porous [.5] burners, volumetric solar receivers, and fibrous and foam insulations; these materials exhibit combined conduction and radiation problems, in which large changes in temperature and thus large variations in thermal conductivity occur [16, 17]. Variable thermal conductivity in convective heat transfer problems may also be noticed in the heat exchangers and cooling systems of electronic devices. In the simulation of these cases, assuming the constancy of thermal conductivity may result in considerable error.

        In this work, the LBM was employed for the first time to analyze non-linear, non-Fourier, two-dimensional (2D) heat transfer with variable thermal conductivity. After the convergence and accuracy of the method was ascertained, the temperature distributions in a 2D plate were obtained under constant thermal conductivity and heat generation conditions. The derived temperature distributions were then compared.

2. Formulation

        The 2D geometry under consideration and its boundary conditions are illustrated in Fig. 1. Four boundaries of the plate were kept at constant temperature, and the wall at x=0.5 m was assumed to be adiabatic. The equation for non-Fourier conduction when two space dimensions are considered in analysis can be written as

           Conduction Heat transfer in a 2-D geometry can be described as Eq. (4)

        Substituting q from Eq. (3) in Eq. (4) yields

where , , k, and g are the density, specific heat coefficient, thermal conductivity, and volumetric heat generation, respectively. In Eq. (5), if thermal lag (relaxation time) , it takes a form governed by Fourier heat conduction. Applying the LBM for simulation necessitates the use of non-dimensional equations. To non-dimensionalize a system of equations, the following steps must be implemented:

1. Identify all independent and dependent variables [.6] 
2. Replace each of the variables with a quantity scaled relative to a characteristic unit of measure to be determined
3. Divide the obtained equation in step 2 by the [.7] coefficient of the highest-order polynomial or derivative term.
4. Judiciously choose a definition of the characteristic unit for each variable so that the coefficients of as many terms as possible become 1.

For generalization, non-dimensional time , distances  and , temperature , heat flux , and volumetric heat generation  are defined in Eq. (6) in the manner recommended in [18]:

[.8] 

        The non-dimensional forms of Eqs. (3) and (4) are

      One of the most common approaches to considering the effects of source terms is attaching an additional term to the standard LBM, as was done in [19]. This approach is reflected in the equation below:

where  is the particle velocity distribution function,  is the weight function, S denotes the source term,  represents the relaxation factor,  is the discrete particle vector,  is the lattice grid, and  is the time increment. Eq. (9) is very convenient to use and accurate for steady source term situations [20].

        For the geometry under consideration, the  (Fig. 2) lattice was used. In this lattice, velocity  and the corresponding weight function  are given as Eqs. (10) and (11):

Fig. 2. Nine-particle velocity in the 2D geometry lattice

        The solution of Eq. (10) in terms of  yields non-dimensional temperature distribution  in the plate as follows:

        Using non-dimensional temperature distribution , non-dimensional heat flux  can be calculated by solving Eq. (3). If the LBM is employed,  is simply obtained according to Eq. (13):

where  and  are the representative horizontal and vertical axes, respectively (Fig. 1). The solution of Eq. (9) needs an equilibrium distribution function , which, for this problem, is assumed to take the following form:

        The equilibrium distribution function must satisfy Eqs. (12) and (13). The satisfaction of the equations is expressed thus:

        Given that Eqs. (8) and (9) are equivalent and represent governing equations of the same problem, they should provide the same results under macroscopic and mesoscopic approaches. The classical Chapman–Enskog expansion is handled to confirm that the results of the mesoscopic approach accord with those of the macroscopic method. [.9] Accordingly, the source term (S) in Eq. (9) can be achieved through this expansion analysis, as reflected in Eq. (16). Details on obtaining source terms are provided in Appendix A.

        Chapman–Enskog multi-scale expansion can be used to confirm that the relationship between relaxation time  and thermal diffusivity  is obtained via Eq. (17):

        In this study, the thermal conductivity coefficient was considered the quadratic function of non-dimensional temperature, as demonstrated in Eq. (18),

2. 1. Boundary Conditions

        The boundaries of the 2D plate were positioned on the lateral sides of the computational domain. In the upper, right, left, and bottom boundaries, , , , and respectively, are unknown and should be specified for streaming[.10] . The boundary conditions are shown in Fig. 1. The unknown distribution functions can be regarded as Eqs. (19) to (24):

To solve the governing equation, (i.e., Eq. (9)), in addition to, the boundary conditions, the initial condition requires identification. This condition is written as follows:

3. Results and Discussion

        Grid convergence reflects result improvement derived through effective small cell sizes for calculations[.11] . A calculation should come on the correct answer as a mesh becomes finer, hence the term “grid convergence”. To do this, grid-independent numerical results for the examined 2D plane should be obtained using several different grid resolutions. As shown in Fig. 3, the number of grids adopted in this work ranged from 5850 to 31850. With 23400 grids, the temperature difference between two consecutive time periods did not exceed .

        In any practical application, a number of important decisions need to be made. These include decisions regarding the number of iterations, the spacing between the iterations retained for the final analysis, and the number of initial burns in discarded iterations. In the problem pursued in this work, the convergence of results was achieved by increasing the number of iterations in the LBM from 1000 to 16000. Fig. 4 presents the temperature distribution at y=0.4 m along the x direction. After 15000 iterations, the results coincided. 

Fig. 5. Comparison of centerline  temperature distributions in the present study and [23].

To validate the results, we analyzed non-Fourier heat conduction in the 2D plate under known temperatures in the four boundaries and compared the temperature distributions at steady-state conditions in the centerline  with those reported in the literature (i.e., [23]). Fig. 5 indicates excellent agreement between the findings.  

        The thermal perturbation source is centrally located on the north and west boundaries of the plate and influences the temperature distribution along the plate[.12] . Figs. 6 and 7 show temperature distributions at x=0.25 m, 0.5 m, and 0.7 m at dimensionless time periods [.13]  respectively. Over time, temperature along the considered lines increased to reach a steady state.

        Temperature distributions in the plate under heat generation conditions at dimensionless time periods , respectively, and the steady state for different locations are shown in Fig. 6. The figure clearly shows the hyperbolic nature of the problem. At , only the limited regions close to the boundaries were affected. With the passage of time, the influence propagated in all parts of the [.14] 2D geometry domain[.15] . At  or when steady-state conditions were reached, the hyperbolic behavior died out, and the conditions became identical to those observed in Fourier conduction.

        Figs. 8 and 9 illustrate the effects of volumetric heat generation on different parts of the plate. In this case, a heat source was assumed to exist at one part of the plate, and other boundaries were assumed to be located at specified boundary conditions. The non-dimensional volumetric heat generation was set at . Heat generation exerted a minimal effect in the beginning compared with the effect imposed under steady-state conditions. This result is attributed to the fact that heat generation entails some time to influence the temperature profile[.16] . Temperature distribution under constant thermal conductivity was reviewed (Fig. 8).A difference was found between the heat generation conditions and the conditions wherein thermal conductivity was a function of temperature. This finding is due to the fact that temperature increased, according to Eq. (27), as thermal conductivity increased. Given the material’s considerable effectiveness in transferring heat energy under heat generation conditions, the temperature gradient of the plate decreased to a level lower than that observed under constant thermal conductivity.

4. Conclusion [.17] 

The LBM was applied in this work to solve the energy equation of a problem with temperature-dependent thermal conductivity. This study is the first to use the LBM in solving the non-Fourier heat conduction problem in a 2D geometry with uniform lattices. The source terms in the LBM formulation were obtained using Chapman–Enskog expansion. Temperature distributions at four temporal periods, including the steady state when the plate was subjected to heat generation and constant thermal conductivity conditions, were also examined. The results showed that because of changes in the ability of the material to transfer heat energy, the temperature gradients under constant and variable thermal conductivity differed.