Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 55, 1, pp. 167-175, Warsaw 2017 DOI: 10.15632/jtam-pl.55.1.167 MODELING OF THERMAL PROCESSES PROCEEDING IN A 1D DOMAIN OF CRYSTALLINE SOLIDS USING THE LATTICE BOLTZMANN METHOD WITH AN INTERVAL SOURCE FUNCTION Alicja Piasecka-Belkhayat, Anna Korczak Silesian University of Technology, Institute of Computational Mechanics and Engineering, Gliwice, Poland e-mail: alicja.piasecka@polsl.pl; anna.korczak@polsl.pl The interval lattice Boltzmann method (ILBM) with an uncertainly defined internal heat source function is used to simulate heat transfer in a thin silicone film. The solution to the interval Boltzmann transport equations has been obtained taking into account the rules of directed interval arithmetics. A similar analysis has been done using the sensitivity model where the Boltzmann transport equations and boundary-initial conditions have been diffe- rentiatedwith respect to the no-interval heat source value. The knowledge of the sensitivity function distribution and the application of the Taylor formula allow one to find the border solutions of the problem analyzed, which (to some extent) correspond to the solution obta- ined under the assumption of the uncertainly defined source function. In the final part of the paper, numerical computations obtained for both methods are presented. Keywords: lattice Boltzmannmethod, directed interval arithmetics, sensitivity analysis 1. Introduction The problem of heat transfer in nano-layers is frequently encountered in many fields of science and engineering such as mechanical engineering, thermal management of electronic cooling and improvement of performance of heat transfer systems (Escobar et al., 2006; Huanga et al., 2005; Joshi andMajumdar, 1993;Mansoor andYilbas, 2011, 2014).Heat transfer problems are usually solved using equations with deterministic thermophysical parameters (Eshraghi and Felicelli, 2012; Narumanchi et al., 2003). However, in most cases of the engineering practice, values of these parameters cannot be defined with a high precision and, in such cases, it is much more convenient to define these parameters as interval numbers (Piasecka-Belkhayat and Korczak, 2014, 2016). In this paper, an interval version of the lattice Boltzmann method with the uncertainly defined heat source function has been presented with the application of the directed interval arithmetics.Thesolutionobtainedcorrespondsto±5%perturbationsof theheat source function. The results of numerical computations (energy and temperature heating curves at the selected points) have an interval form, of course. Additionally, the sensitivity analysis with respect to the constant heat source function has been done (Chonga et al., 2016; Dems and Rousselet, 1999; Goethals et al., 2011; Hwang et al., 2016). Theheat source value has been assumedas themiddle value of the heat source interval. The application of the sensitivity function distribution and the Taylor formula with an increment of the source function equal to the half of the width of the heat source interval allows one to find the solution to the boundary-initial problem similar to the solution with some “uncertainties” appearing in the mathematical model. The aim of the paper is comparison of the results obtained using both methods. 168 A. Piasecka-Belkhayat, A. Korczak 2. Directed interval arithmetics Let us consider a directed interval awhich can be defined as a set D of all directed pairs of real numbers defined as follows (Neumaier, 1990; Piasecka-Belkhayat, 2011a,b) a= [a−,a+] for a−,a+ ∈ R (2.1) where a− and a+ denote the beginning and the end of the interval, respectively. The left or the right endpoint of the interval a can be denoted as as, s ∈ {+,−}, where s is a binary variable. This variable can be expressed as a product of two binary variables and is defined as ++=−−=+ +−=−+=− (2.2) An interval is called proper if a− < a+, improper if a− > a+ and degenerate if a− = a+. The set of all directed interval numbers can be written as D = P∪ I, where P denotes a set of all directed proper intervals and I denotes a set of all improper intervals. Additionally, a subset Z = ZP ∪ZI ⊂ D should be defined, where ZP = {a∈ P| a − ¬ 0¬ a+} ZI = {a∈ I| a + ¬ 0¬ a−} (2.3) For directed interval numbers, two binary variables are defined. The first of them is the direction variable τ(a)= { + if a− ¬ a+ − if a− >a+ (2.4) and the other is the sign variable σ(a)= { + if a− > 0, a+ > 0 − if a− < 0, a+ < 0 a∈ D\Z (2.5) For the zero argument σ([0,0]) = σ(0) =+, for all intervals including the zero element a∈ Z, σ(a) is not defined. The sum of two directed intervals a= [a−,a+] and b= [b−,b+] can be written as a+ b= [a−+ b−,a++ b+] a,b∈ D (2.6) The difference is of the form a− b= [a−− b+,a+− b−] a,b∈ D (2.7) The product of the directed intervals is described by the formula ab=                          [a−σ(b)b−σ(a),aσ(b)bσ(a)] a,b∈ D\Z [aσ(a)τ(b)b−σ(a),aσ(a)τ(b)bσ(a)] a∈ D\Z, b∈ Z [a−σ(b)bσ(b)τ(a),aσ(b)bσ(b)τ(a)] a∈ Z, b∈ D\Z [min(a−b+,a+b−),max(a−b−,a+b+)] a,b∈ ZP [max(a−b−,a+b+),min(a−b+,a+b−)] a,b∈ ZI 0 (a∈ ZP,b∈ ZI) ∨ (a∈ ZI,b∈ ZP) (2.8) Modeling of thermal processes proceeding in a 1D domain... 169 The quotient of two directed intervals can be written as a/b=    [a−σ(b)/bσ(a),aσ(b)/b−σ(a)] a,b∈ D\Z [a−σ(b)/b−σ(b)τ(a),aσ(b)/b−σ(b)τ(a)] a∈ Z, b∈ D\Z (2.9) In the directed interval arithmetics, two extra operators are defined, the inversion of summation −Da= [−a −,−a+] a∈ D (2.10) and the inversion of multiplication 1/Da= [1/a −,1/a+] a∈ D\Z (2.11) So, two additional mathematical operations can be defined as follows a−D b= [a −− b−,a+− b+] a,b∈ D (2.12) and a/Db=    [a−σ(b)/b−σ(a),aσ(b)/bσ(a)] a,b∈ D\Z [a−σ(b)/bσ(b),aσ(b)/bσ(b)] a∈ Z, b∈ D\Z (2.13) Now, it is possible to obtain the number zero by subtraction of two identical intervals a−Da=0 and the number one as the result of division a/Da = 1, which is impossible when applying classical interval arithmetics (Markov, 1995). 3. Boltzmann transport equation The Boltzmann transport equation (BTE) is one of the fundamental equations of solid state physics and takes the following form (Escobar et al., 2006) ∂f ∂t +v∇f = f0−f τr +gef (3.1) where f is the phonon distribution function, f0 is the equilibriumdistribution function given by the Bose-Einstein statistics, v is the phonon group velocity, τr is the relaxation time and gef is the phonon generation rate due to electron-phonon scattering. In order to take advantage of the simplifying assumption of the Debye model, the BTE can be transformed into an equation of the carrier energy density, and for a one-dimensional problem has the following form (Escobar et al., 2006) ∂e ∂t +v∇e=− e−e0 τr +qv (3.2) where e is the phonon energy density, e0 is the equilibrium phonon energy density and qv is the internal heat source related to a unit of volume. Equation (3.2) must be supplemented by the adequate boundary-initial conditions. Using the Debye model, the relation between the phonon energy density and lattice tempe- rature is given by the following formula (Escobar et al., 2006) e(T)= ( 9ηkb Θ3 D ΘD/T ∫ 0 z3 exp(z)−1 dz ) T4 (3.3) 170 A. Piasecka-Belkhayat, A. Korczak whereΘD is the Debye temperature of the solid, kb is the Boltzmann constant, T is the lattice temperature while η is the number density of phonons, and can be calculated using the formula (Escobar et al., 2006) η= 1 6π2 (kbΘD ~ω )3 (3.4) where ~ is the Planck constant divided by 2π and ω is the speed of sound in the analysed material. 4. Interval lattice Boltzmann method The lattice Boltzmannmethod (LBM) is a numerical technique for simulation of fluid flows and heat transfer. Here LBMhas been successfully applied to simulate heat transfer in nano layers. Unlike the conventional numerical methods based on discretizations of macroscopic continuum equations, the LBM is based on nanoscale models and heat transfer equations. In this paper, it is shown how the LBM solves a discretized set of the Boltzmann transport equation (BTE) in the case of interval values appearing in the mathematical model. Then the interval Boltzmann transport equation for a one-dimensional problem has the following form (Piasecka-Belkhayat and Korczak, 2014) ∂e(x,t) ∂t +v ∂e(x,t) ∂x =− e(x,t)−e0(x,t) τr + qv(x,t) (4.1) where e(x,t) is the interval phonon energy density, e0(x,t) is the interval equilibrium phonon energy density, τr is the relaxation time, v1 = v and v2 = −v (see Fig. 1) and qv(x,t) is the interval heat source, x is the spatial coordinate and t is the time. Fig. 1. Directions of the lattice vibrations The interval total phonon energy density is defined as the sumof phonon energy densities in all directions. In the paper, a one-dimensionalmodelwith twodirections of the phononvelocities is assumed (Piasecka-Belkhayat and Korczak, 2014, 2016) e(x,t)= e1(x,t)+e2(x,t)= 2 ∑ d=1 ed(x,t) (4.2) where e1(x,t) is thephonon energydensity in thepositivexdirectionwhile e2(x,t) is the phonon energy density in the negative x direction, dmeans the lattice direction (see Fig. 1). In the interval latticeBoltzmannmethod it isneeded to solve systemof twopartial differential equations allowing one to compute phonon energy in different lattice nodes according to the following equations (Piasecka-Belkhayat andKorczak, 2014) ∂ed(x,t) ∂t +(−1)d−1v ∂ed(x,t) ∂x =− ed(x,t)−e 0 d(x,t) τr +qv(x,t) d=1,2 (4.3) where v=∆x/∆t is the component of velocity along the x-axis,∆x is the lattice distance from site to site, ∆t= tf+1− tf is the time step needed for a phonon to travel from one lattice site to the neighboring lattice site and e0d(x,t)= e(x,t)/d (4.4) Modeling of thermal processes proceeding in a 1D domain... 171 The set of equations (4.3) must be supplemented by the boundary-initial conditions (Goethals et al., 2011) x=0 : e1(0, t) = e(Tb1) x=L : e2(L,t)= e(Tb2) t=0 : e(x,0)= e(T0) (4.5) where Tb1 and Tb2 are the boundary temperatures, T0 is the initial temperature. The approximate form of equations (4.3) is of the following form (e1) f+1 i+1 = ( 1− ∆t τr ) (e1) f i + ∆t τr (e01) f i +∆tqv (e2) f+1 i−1 = ( 1− ∆t τr ) (e2) f i + ∆t τr (e02) f i +∆tqv (4.6) After subsequent computations, the interval lattice temperature is determined according to the rules of directed interval arithmetics using the formula (see Eq. (3.3)) T f = 4 √ √ √ √ √ √efΘ3 D ( 9ηkb ΘD/T f−1 ∫ 0 z3 exp(z)−1 dz ) −1 (4.7) 5. Sensitivity analysis In order to analyze the sensitivity of the phonon energy density field, the governing equations should be differentiated with respect to the chosen parameter (Kleiber, 1997). In the paper, the sensitivity analysis is presented with respect to the value of the internal heat source. The Boltzmann transport equation for the one-dimensional problem and the constant value of the heat source qv has the following form (Kałuża et al., 2016; Majchrzak and Mochnacki, 2014; Mochnacki andMajchrzak, 2007; Mohebbi and Sellier, 2016) ∂ed(x,t) ∂t +(−1)d−1v ∂ed(x,t) ∂x =− ed(x,t)−e 0 d(x,t) τr + qv d=1,2 (5.1) with the boundary-initial conditions x=0 : e1(0, t) = e(Tb1) x=L : e2(L,t)= e(Tb2) t=0 : e(x,0)= e(T0) (5.2) Using the direct approach of sensitivity analysis, equation (5.1) is differentiated with respect to qv (Jasiński, 2014; Mochnacki andMajchrzak, 2007; Mohebbi and Sellier, 2016) ∂ ∂qv (∂ed(x,t) ∂t ) +(−1)d−1v ∂ ∂qv (∂ed(x,t) ∂x ) =− 1 2τr ∂ed(x,t) ∂qv + ∂qv ∂qv d=1,2 (5.3) Next, differentiation of boundary-initial conditions (5.2) leads to the following formulas x=0 : ∂e1(0, t) ∂qv = ∂e(Tb1) ∂qv =0 x=L : ∂e2(L,t) ∂qv = ∂e(Tb2) ∂qv =0 t=0 : ∂e(x,0) ∂qv = ∂e(T0) ∂qv =0 (5.4) 172 A. Piasecka-Belkhayat, A. Korczak To equation (5.3) and boundary-initial conditions (5.4), the sensitivity functions Ud(x,t,qv) = ∂ed(x,t)/∂qv are introduced ∂Ud(x,t,qv) ∂t +(−1)d−1v ∂Ud(x,t,qv) ∂x =− 1 2τr Ud(x,t,qv)+1 d=1,2 x=0 : U1(0, t,qv)= 0 x=L : U2(L,t,qv)= 0 t=0 : U(x,0,qv)= 0 (5.5) Then equations (5.3) are the following ∂Ud(x,t,qv) ∂t +(−1)d−1v ∂Ud(x,t,qv) ∂x =− Ud(x,t,qv)−U 0 d(x,t,qv) τr +1 d=1,2 (5.6) where U0d(x,t,qv)= ∂e0d(x,t) ∂qv = ∂ ∂qv (e(x,t) 2 ) = U(x,t,qv) 2 (5.7) while U(x,t,qv)= 2 ∑ d=1 Ud(x,t,qv) (5.8) The phonon energy density function e(x,t,qv ±∆qv) is expanded into the Taylor series taking into account the first two components according to the formulas e(x,t,qv +∆qv)≈ eb(x,t)+ ∂e(x,t) ∂qv ∆qv e(x,t,qv −∆qv)≈ eb(x,t)− ∂e(x,t) ∂qv ∆qv (5.9) where∆qvis a certain incrementof the source function, andthe startingpointeb(x,t) corresponds to the basic solution. Taking into account the sensitivity functionU(x,t,qv)= ∂e(x,t)/∂qv, one obtains e(x,t,qv +∆qv)≈ eb(x,t)+U(x,t,qv)∆qv e(x,t,qv −∆qv)≈ eb(x,t)−U(x,t,qv)∆qv (5.10) and a certain increment of the energy function∆e can be calculated using the formula ∆e(x,t)≈ 2U(x,t,qv)∆qv (5.11) 6. Numerical example In the paper, heat transfer in a one-dimensional silicon film of dimensionL=200nm has been analyzed. The following input data have been introduced: relaxation time τr = 6.53ps, Debye temperature ΘD = 640K, initial temperature T0 = 300K, boundary conditions Tb1 = Tb2 = 300K, lattice distance∆x=20nm and the time step∆t=5ps. In the first example, the interval value of the heat source function has been considered qv = [10 18−0.05 ·1018,1018+0.05 ·1018]W/m3. Figure 2a illustrates the interval heating curves of the phonon energy at the internal nodes 1 (20nm), 2 (80nm) and 3 (140nm). Modeling of thermal processes proceeding in a 1D domain... 173 Fig. 2. Energy heating curves: (a) – first method, (b) – secondmethod In the second example, the no-interval value of the heat source function has been introduced qv = 10 18W/m3 and the sensitivity analysis with respect to the heat source parameter has been applied. In this model, an increment of the heat source parameter has been introduced as ∆qv =0.05 ·10 18W/m3. In Fig. 2b, the courses of heating curves of the phonon energy at the same internal nodes are presented. One can see that the both results are similar. Figures 3a and3b present the courses of heating curves taking into account the same internal nodes for the first and second example, respectively. Fig. 3. Temperature heating curves: (a) – first method, (b) – secondmethod 7. Conclusions In the paper, heat transfer in one-dimensional crystalline solids is considered. Themain subject of the paper is the comparison of the results obtained using two methods. In the first method, the interval lattice Boltzmannmethodwith an uncertainly defined internal heat source function is used. The solution to the interval Boltzmann transport equations has been obtained taking into account the rules of directed interval arithmetics. In the second method, the sensitivity 174 A. Piasecka-Belkhayat, A. Korczak analysis with respect to the internal heat source parameter has been done. The application of the sensitivity functions and the Taylor formula enables one to find a solution similar to the solution received using the interval lattice Boltzmann method. References 1. Chonga W.T., Al-MamoonaA., Poha S.Ch., Sawb L.H., Shamshirbandc S.,Mojumder J.Ch., 2016, Sensitivity analysis of heat transfer rate for smart roof design by adaptive neuro-fuzzy technique,Energy and Buildings, 124, 112-119 2. Dems K., Rousselet B., 1999, Sensitivity analysis for transient heat conduction in a solid body – Part I, Structural Optimization, 17, 36-45 3. Escobar R.A., Ghai S.S., Jhon M.S., Amon C.H., 2006,Multi-length and time scale thermal transport using the lattice Boltzmann method with application to electronics cooling, Journal of Heat and Mass Transfer, 49, 97-107 4. Eshraghi M., Felicelli S.D., 2012, An implicit lattice Boltzmann model for heat conduction with phase change, International Journal of Heat and Mass Transfer, 55, 2420-2428 5. Goethals K., Breeschb H., Janssens A., 2011, Sensitivity analysis of predicted night cooling performance to internal convective heat transfer modeling,Energy and Buildings, 43, 2429-2441 6. Huanga S.M., Sun Z., Luk’yanchuk B.S., Hong M.H., Shi L.P., 2005, Nanobump arrays fabricated by laser irradiation of polystyrene particle layers on silicon,Applied Physics Letters, 86, 161911 7. Hwang S., Son Ch., Seo D., Rhee D.H., Cha B., 2016, Comparative study on steady and unsteady conjugate heat transfer analysis of a high pressure turbine blade,Applied Thermal Engi- neering, 99, 765-775 8. Jasiński M., 2014,Modeling of tissue thermal injury formation process with application of direct sensitivity method, Journal of Theoretical and Applied Mechanics, 52, 4, 947-957 9. Joshi A.A.,Majumdar A., 1993, Transient ballistic and diffusive phonon heat transport in thin films, Journal of Applied Physics, 74, 1, 31-39 10. Kałuża G., Majchrzak E., Turchan Ł., 2016, 1D generalized dual-phase lag equation. Sen- sitivity analysis with respect to the porosity, Journal of Applied Mathematics and Computational Mechanics, 15, 1, 49-58 11. Kleiber M., 1997,Parameter Sensitivity in Non-linear Mechanics, J.Willey & Sons, London 12. MajchrzakE.,MochnackiB., 2014,Sensitivity analysis of transient temperaturefield inmicro- domainswith respect to the dual phase lagmodel parameters, International Journal for Multiscale Computational Engineering, 12, 1, 65-77 13. Mansoor S. Bin, Yilbas B.S., 2011, Laser short-pulse heating of silicon-aluminum thin films, Optical and Quantum Electronics, 42, 601-618 14. Mansoor S. Bin, Yilbas B.S., 2014, Phonon transport in aluminum and silicon film pair: laser short-pulse irradiation at aluminum film surface,Canadian Journal of Physics, 92, 12, 1614-1622 15. Markov S.M., 1995, On directed interval arithmetic and its applications, Journal of Universal Computer Science, 1, 514-526 16. Mochnacki B., Majchrzak E., 2007, Identification of macro and micro parameters in solidifi- cationmodel,Bulletin of the Polish Academy of Sciences, Technical Sciences, 55, 1, 107-113 17. Mohebbi F., SellierM., 2016,Estimation of thermal conductivity, heat transfer coefficient, and heat flux using a three dimensional inverse analysis, International Journal of Thermal Sciences, 99, 258-270 Modeling of thermal processes proceeding in a 1D domain... 175 18. Narumanchi S., Murthy J.Y., Amon C.H., 2003, Simulation of unsteady small heat source effects in sub-micron heat conduction, Journal of Heat Transfer, 123, 896-903 19. NeumaierA., 1990, IntervalMethods for System of Equations, CambridgeUniversityPress,Cam- bridge, NewYork, Port Chester, Melbourne, Sydney 20. Piasecka-Belkhayat A., 2011, Interval boundary element method for 2D transient diffusion problem using directed interval arithmetic, Engineering Analysis with Boundary Elements, 35, 3, 259-263 21. Piasecka-BelkhayatA., 2011, Intervalboundaryelementmethod for transientdiffusionproblem in two-layered domain, Journal of Theoretical and Applied Mechanics, 49, 1, 265-276 22. Piasecka-Belkhayat A., Korczak A., 2014, Modelling of transient heat transport in one- -dimensional crystalline solids using the interval lattice Boltzmann method, Recent Advances in Computational Mechanics, T. Łodygowski, J. Rakowski and P. Litewka (Eds.), Taylor & Francis Group, A Balkema Book, London, 363-368 23. Piasecka-Belkhayat A., Korczak A., 2016, Numerical modelling of the transient heat trans- port in a thin gold film using the fuzzy lattice Boltzmannmethod with α-cuts, Journal of Applied Mathematics and Computational Mechanics, 15, 1, 123-135 Manuscript received May 16, 2016; accepted for print July 4, 2016