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

Document Type : Full Lenght Research Article


University of Kashan


The lattice Boltzmann method (LBM) was used to analyze two-dimensional (2D)
non-Fourier heat conduction with temperature-dependent thermal conductivity. To this end, the evolution of wave-like temperature distributions in a 2D plate was obtained. The temperature distributions along certain parts of the plate, which was subjected to heat generation and constant thermal conductivity conditions, were also derived and compared. The LBM results are in good agreement with those reported in other works. Additionally, the temperature contours at four different times in which steady state conditions can be achieved were analyzed. The results showed that thermal conductivity increased with rising temperature. 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.

Keywords: Non-Fourier conduction, lattice Boltzmann method, variable thermal conductivity, constant thermal conductivity, heat generation


Main Subjects

  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]:





        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. 1. Geometry and boundary conditions of the problem under consideration.




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. 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 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.









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









Fig. 7. Contours of temperature distribution derived via the current LBM model at different time periods: (a) , (b) , (c) , and (d)  (steady state).









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. Temperature-dependent 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] 




microscopic velocity vector


speed of thermal wave


specific heat  


propagation velocity in direction i in the lattice


particle distribution function in the i direction


thermal conductivity  


reference length (m)




conduction heat flux


position vector


source term




time (s)

X, Y

dimensions of the 2D rectangular enclosure

Greek symbols


thermal diffusivity coefficient 


non-dimensional temperature


non-dimensional time


relaxation time


direction cosine


time lag


non-dimensional heat flux


i, j

coordinates of a cell center



x, y

for x-y faces of control volume


initial value



Non-dimensional quantities



[1] J. Zhou, Y. Zhang, J.K. Chen, Non-Fourier heat conduction effect on laser induced thermal damage in biological tissues,  Heat Transfer, A 54 (1), 65-70, (2008).        
[2] J.Y. Lin, The non-Fourier effect on the fin performance under periodic thermal conditions, Appl. Math, Model. 22, 629–640, (1998).
[3] D.W. Tang, N. Araki, Non-Fourier heat conduction behavior in finite mediums under pulse surface heating, Mater. Sci. Eng, A 292, 173–178, (2000).
[4] J.R. Ho, C.P. Kuo, W.S. Jiaung, C.J. Twu, Lattice Boltzmann scheme for hyperbolic heat conduction, Numer. Heat Transfer, B 41, 591–607, (2002).
[5] S.C. Mishra, A. Lankadasu, K. Beronov, Application of the lattice Boltzmann method for solving the energy equation of a 2-D transient conduction radiation problem, Int. J. Heat Mass Transfer, 48, 3648–3659, (2005).
[6] S.C. Mishra, H.K. Roy, Solving transient conduction-radiation problems using the lattice Boltzmann method and the finite volume method, J. Compute. Phys. 233, 89–107, (2007).
[7] M. H. Rahimian, I. Rahbari, F. Mortazavi, High order numerical simulation of non-Fourier heat conduction: An application of numerical Laplace transform inversion, International Journal of Heat and Mass Transfer, vol. 51, 51–58, (2014).
[8] W. Dreyer, S. Qamar, Kinetic flux-vector splitting schemes for the hyperbolic heat conduction, J. Comput. Phys. 198 (2), (2004).
[9] H. Chen, J. Lin, Numerical analysis for hyperbolic heat conduction, Int. J. Heat Mass Transfer 36 (11), 2891–2898, (1993).
[10] J.I. Frankel, B. Vick, M.N. Özisik, General formulation and analysis of hyperbolic heat conduction in composite media, Int. J. Heat Mass Transfer 30, 1293–1305, (1987).
[11] B. Abdel-Hamid, Modelling non-Fourier heat conduction with periodic thermal oscillation using the finite integral transform, Appl. Math, Model 23, 899–914, (1999).
[12] T.M. Chen, C.C. Chen, Numerical solution for the hyperbolic heat conduction problems in the radial–spherical coordinate system using a hybrid Green’s function method, Int. J. Therm. Sci. 49, 1193–1196, (2010).
[13] X. Lu, P. Tervola, M. Viljanen, Transient analytical solution to heat conduction in composite circular cylinder, Int. J. Heat Mass Transfer, 49, 341–348, (2006).
[14] G.E. Cossali, Periodic heat conduction in a solid homogeneous finite cylinder, Int. J. Therm. Sci. 48, 722–732, (2009).
[15] A. Moosaie, Axsymmetric non-Fourier temperature field in a hollow sphere, Arch. Appl. Mech. 79, 679–694, (2009).
[16] H. Ahmadikia and M. Rismanian, Analytical solution of non-Fourier heat conduction problem on a fin under periodic boundary conditions, Journal of Mechanical Science and Technology, vol. 25 (11) 2919-2926, (2011).
[17] R. Siegel, J. Howell, Thermal Radiation heat Transfer, Taylor & Francis: New York, (2002).
[18] C. Mishra, H. Sahai, Analyses of non-Fourier heat conduction in 1-D cylindrical ans spherical geometry- An application of the lattice Boltzmann method, International Journal of Heat and Mass Transfer, 55, 7015-7023, (2012).
[19] Q. Zou, S. Hou, G. D. Doolen, Analytical solutions of the lattice Boltzmann BGK model, Journal of Statistical Physics, vol(81), 319–334, (1995).
[20] P. Lallemand, L. S. Luo, heory of the lattice Boltzmann method-acoustic and thermal properties in two and three dimensions, Physical Review, E68, (036706), 1–25, (2003).
[21] T. M. Chen, C. C. Chen, Numerical solution for the hyperbolic heat conduction problems in the radial-spherical coordinate system using a hybrid Green’s function method, International Journal of Thermal Sciences, (2010).
[22] C. C. Wang, Direct and inverse solutions with non-Fourier effect on the irregular shape, International Journal of Heat and Mass Transfer, vol (53), (13-14), 2685–2693, (2010).
[23] C. Mishra, B. Mondal, T. Kush, B. Sima Rama Krishna, Solving transient heat conduction problems on uniform and non-uniform lattices using the lattice Boltzmann method, International Communications in Heat and Mass Transfer, vol(36), 322–328, (2009).
  • Receive Date: 26 October 2016
  • Revise Date: 04 September 2017
  • Accept Date: 04 October 2017
  • First Publish Date: 01 May 2018