Document Type : Full Lenght Research Article
Authors
^{1} University of Kashan
^{2} University of Kashan
Abstract
Keywords
Main Subjects
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 nonFourier effect is a more reliable means through which to explain the aforementioned [.1] processes. NonFourier motion (visually perceived motion that cannot be explained simply on the basis of the auto correlational structure of a visual stimulus) is a wellrecognized phenomenon that is generally considered to be due to the nonlinear 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):
, 
(1) 
Most thermal conduction or mass diffusion problems are depicted and investigated using Fourier’s law, on whose basis nonphysical 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 nonFourier heat conduction theories. The mathematical statement of nonFourier 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] :
(2) 
where is the thermal lag time, that is, the rate of thermal diffusivity ( ) to the square of the speed of heat wave propagation , .
The nonFourier 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 nonFourier heat conduction mechanism, there are few studies about it. Ho et al. [4] and Mishra et al. [5, 6], investigated nonFourier heat conduction in onedimensional planar geometry. In some studies also analyzed nonFourier 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 nonlinear, nonFourier, twodimensional (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 nonFourier conduction when two space dimensions are considered in analysis can be written as
(3) 
Conduction Heat transfer in a 2D geometry can be described as Eq. (4)
(4) 
Substituting q from Eq. (3) in Eq. (4) yields
(5) 
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 nondimensional equations. To nondimensionalize a system of equations, the following steps must be implemented:
For generalization, nondimensional time , distances and , temperature , heat flux , and volumetric heat generation are defined in Eq. (6) in the manner recommended in [18]:

(6) 
The nondimensional forms of Eqs. (3) and (4) are

(7a) 
(7b) 


(8) 
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:
(9) 
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):

(10) 
Fig. 1. Geometry and boundary conditions of the problem under consideration.

Fig. 2. Nineparticle velocity in the 2D geometry lattice

(11) 
The solution of Eq. (10) in terms of yields nondimensional temperature distribution in the plate as follows:
(12) 
Using nondimensional temperature distribution , nondimensional heat flux can be calculated by solving Eq. (3). If the LBM is employed, is simply obtained according to Eq. (13):
(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:

(14) 
The equilibrium distribution function must satisfy Eqs. (12) and (13). The satisfaction of the equations is expressed thus:
(15) 
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.
(16) 
Chapman–Enskog multiscale expansion can be used to confirm that the relationship between relaxation time and thermal diffusivity is obtained via Eq. (17):

(17) 
In this study, the thermal conductivity coefficient was considered the quadratic function of nondimensional temperature, as demonstrated in Eq. (18),
(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):
(19) 
(20) 
(21) 
(22) 
(23) 
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:

(24) 
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, gridindependent 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. 3. Centerline (y=0.4 m) temperature evolution at different numbers of grids. 
Fig. 4. Centerline (y=0.4 m) temperature evolution at different numbers of iterations. 
Fig. 5. Comparison of centerline temperature distributions in the present study and [23].
To validate the results, we analyzed nonFourier heat conduction in the 2D plate under known temperatures in the four boundaries and compared the temperature distributions at steadystate 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 steadystate 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 nondimensional volumetric heat generation was set at . Heat generation exerted a minimal effect in the beginning compared with the effect imposed under steadystate 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 temperaturedependent thermal conductivity. This study is the first to use the LBM in solving the nonFourier 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.
(a) 
(b) 

(c) 
Fig. 6. Temperature distributions at different time periods and (a) x=0.25 m, (b) x=0.5 m, and (c) x=0.7 m 
(a) 
(b) 
(c) 
(d) 
Fig. 7. Contours of temperature distribution derived via the current LBM model at different time periods: (a) , (b) , (c) , and (d) (steady state). 
(a) 
(b) 
(c) 
(d) 
Fig. 8. Temperature distributions at different time periods and (a) x=0.25 m, (b) x=0.5 m, (c) x=0.7 m, and (d) y=0.4 m. Temperaturedependent thermal conductivity is denoted by a solid line; constant thermal conductivity is denoted by a solid line with a shaded circle on it; and heat generation is denoted by a dashed line.[M18] 
Nomenclature 

microscopic velocity vector 

C 
speed of thermal wave 
specific heat 

propagation velocity in direction i in the lattice 

particle distribution function in the i direction 

k 
thermal conductivity 
L 
reference length (m) 
m 
direction 
q 
conduction heat flux 
position vector 

S 
source term 
T 
temperature 
t 
time (s) 
X, Y 
dimensions of the 2D rectangular enclosure 
Greek symbols 

thermal diffusivity coefficient 


nondimensional temperature 
nondimensional time 

relaxation time 

direction cosine 

time lag 


nondimensional heat flux 
Subscripts 

i, j 
coordinates of a cell center 
ref 
reference 
x, y 
for xy faces of control volume 
0 
initial value 
Superscripts 

* 
Nondimensional quantities 