Energy transport across submicron porous structures: A Lattice Boltzmann study

International Journal of Heat and Mass Transfer 72 (2014) 479–488

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Energy transport across submicron porous structures: A Lattice Boltzmann study Ankur Chattopadhyay, Arvind Pattamatta ⇑ Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600036, India

a r t i c l e

i n f o

Article history: Received 11 November 2013 Received in revised form 17 January 2014 Accepted 18 January 2014 Available online 4 February 2014 Keywords: Lattice Boltzmann Method (LBM) Boltzmann Transport Equation (BTE) Phonon Electron Semiconductor nanopore Metallic nanopore

a b s t r a c t The technological advancement in the field of energy conversion demand exploring new materials with low value of thermal conductivity. The reduction in the thermal conductivity of a material can be achieved by decreasing its device dimension and increasing the porosity of the material. Lattice Boltzmann Method (LBM), a discrete form of Boltzmann Transport Equation (BTE) has been employed in the present study to examine the energy transport across semiconductor and metal structures. A newly proposed lattice weight system is exercised to alleviate the anisotropic phonon movement in the D2Q9 lattice. The study explores the comparison of different pore configurations classified based on shape and material with varying interfacial pore scattering. In semiconductors, the selection of pore geometry and specularity parameter has been exploited comprehensively in scaling down the thermal properties. Results reveal for rectangular and square pores, thermal conductivity varies inversely with the pore scattering area. Also the comparative analysis of Debye and Sine dispersion model is performed. The role of electron–phonon contribution towards thermal conductivity in metals is analyzed. When energy is transported across the metallic gold film, the electron contribution of the total thermal conductivity exhibits considerable size effects, however phonon contribution remains unaltered with that of its bulk value. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction The evaluation of thermal conductivity of micro/nano porous structures has garnered significant attention in recent years. While determining the performance of the thermo-electric (TE) materials using the figure of merit ZT = S2rT/k (where S is the Seebeck coefficient, r is the electrical conductivity, T is the absolute temperature and k is thermal conductivity), the estimation of the thermal properties of a material becomes a key challenge. The higher value of ZT can be achieved by controlling the thermal conductivity rather than manipulating the electrical conductivity [1]. Therefore, from the energy conversion point of view, the contribution to the low value of thermal conductivity is essential for the development of efficient TE materials. A few studies have already explored the evaluation of thermal conductivity of the porous material using both experimental and computational techniques. Benedetto et al. [2] measured the thermal conductivity of porous silicon layers of various thicknesses belonging to sub continuum scale with varying porosity and wafer type. Gesele et al. [3] also reported the size effects in porous silicon membranes of variable doping level over a wide temperature range. Drost et al. [4] observed the reduction of thermal ⇑ Corresponding author. Tel.: +91 4422574654. E-mail address: [email protected] (A. Pattamatta). 0017-9310/$ - see front matter Ó 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.01.040

conductivity in porous silicon layers due to oxidation. Also other experimental investigations [5,6] reflected that the measured thermal conductivity of porous silicon was considerably lower than that of its bulk value. As far as computations are concerned, Chung et al. [7] incorporated the fundamental Boltzmann Transport Equation (BTE) in analyzing the thermal conductivity of porous silicon and revealed the necessity of considering pore scattering and pore randomness. Besides, the phonon transport in micro-nano porous structures has been examined by other numerical techniques like Molecular Dynamics (MD), Monte Carlo (MC) simulations, etc. [8–10]. In fact Prasher [11,12] also developed specific analytical models in determining thermal conductivity of micro-nano porous structures and composites. However in recent times a few studies [13–17] have experimented with an alternate method called the Lattice Boltzmann Method (LBM) while solving phonon transport. LBM, originated from the discretized version of BTE for modeling energy transport offers reasonable computational simplicity compared to BTE. LBM has been extensively promoted in the study of flow physics at the continuum level [18,19] as well as micro scale [20]. A few studies have dealt with the application of LBM to examine the phonon transport in various applications from electronic cooling to data storage [13,14]. Recently, mode dependent phonon transport has been investigated by Heino [15]. A study where LBM is coupled with a finite difference model to analyze the macro/micro heat

480

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

Nomenclature kB L D(x)

veff

eeq i Q Kn p f e000 feq v  h x T k t e C

j

Boltzmann’s constant (J/K) characteristic length density of states (1/m3) effective phonon group velocity energy equilibrium distribution (J/m3) heat flux (W/m2) knudsen number momentum (kg m/s) phonon distribution function phonon energy density (J/m3) phonon equilibrium distribution function phonon group velocity (m/s) reduced Planck’s constant (J s) spatial variable (m) temperature (K) thermal conductivity (W/m K) time (s) total phonon energy distribution (J/m3) total specific heat (J/m3 K) wave vector

Greek symbols h,/ angles r gradient operator K mean free path (MFP) x phonon frequency (rad/s) s phonon relaxation time (s) Subscripts col collision or scattering term eff effective group velocity el electron eq equilibrium distribution function 2 location of cold temperature boundary 1 location of hot temperature boundary n normal direction ph phonon tot total energy distribution

transfer is also performed by Christensen et al. [16]. Nabovati et al. [17] analyzed two dimensional lattice structures and proposed that D2Q7 yields better performance than D2Q9 without considering the center point of a particular lattice. It is quite interesting to note that most of the recent studies [13–17] have exercised LBM for solving phonon energy transport in high Knudsen number regimes. In spite of these studies, there is no literature available to the best of authors’ knowledge, where the results of submicron energy transport using LBM in porous configurations are evaluated. Therefore in this context, the present study demonstrates the applicability of the LBM solver in predicting the thermal conductivity of porous structures (both semiconductors and metals), especially at submicron regime. A detailed methodology of applying LBM is explained and implemented in numerical examples in the following sections. A novel lattice weight configuration is also proposed in D2Q9 lattice, in order to take care of its inherent non isotropic lattice structure and non uniformity in lattice speeds. Initially, the present study focuses upon validating the present LBM solver with the computational and experimental results available in the literature. Further, this study is extended to compare the Debye and Sine dispersion models for Knudsen numbers ranging between 0.1 and 0.4, in terms of estimating the thermal conductivities with varying parameters like pore configurations, porosity and pore scattering. Moreover the effect of porosity and pore geometries are also examined in metallic gold films where the role of electron and phonon contribution in energy transport is analyzed.

  @f f eq  f ¼ @t col s

The general form of BTE in the absence of any heat generation is expressed as follows [21]:

ð1Þ

where v denotes the phonon group velocity. The distribution function of an ensemble of particles in domain of interest is described by ‘f’ and ‘F’ denotes the external force. The complicated scattering term in the right hand side represents the rate of change of a distribution function due to collision which restores the system to equi-

ð2Þ

where feq represents the equilibrium distribution function that reinstates the system into equilibrium after the collision and s is known as the relaxation time. The equilibrium distribution function can be calculated from the Bose–Einstein relationship:

f eq ¼

exp



1 x h kB T



ð3Þ

1

The phonon energy density can be described from distribution function as:

xfDðxÞ e000 ¼ h

ð4Þ

 is the where e is the phonon energy density per unit solid angle, h reduced Planck’s constant, x is the phonon frequency and D(x) is the density of states. The total equilibrium energy can be calculated by integrating the energy density over the entire frequency band and is related to the product of total heat capacity (C) of the solid and physical temperature (T) as follows: 000

eeq ¼

2. Governing equations

  @f @f @f þ m:rf þ F: ¼ @t @p @t col

librium from an unbalanced state. However in its most stringent form, the scattering term is non-linear and can be simplified by using relaxation time approximation as proposed by Bhatnagar et al. [22]:

1 4p

Z

xfDðxÞ ¼ 1 h 4p Dx

Z Dx

e000 dld/ ¼

CT 4p

ð5Þ

where l is the direction cosine of x axis and / is the angular resolution. In order to reduce the computational cost of solving rigorous BTE, a discretized version of the BTE, named as LBM has also been proposed for solving phonon transport [13–15,17]. The numerical domain is composed of the sequential distribution of lattices. The LBM technique can be classified into two processes – collision and propagation. The phonons travel ballistically with a constant group velocity from one node to successive nodes as guided by lattice directions. There they collide with other phonons generated at the consecutive lattices and continue to March towards the

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

adjacent lattices. The phonons follow similar recurring events at immediate nodes till the attainment of equilibrium. After replacing the distribution function with the energy density defined in Eq. (4), the discrete BTE can be rewritten in the following manner:

ei ðx; t þ DtÞ  ei ðx; tÞ ei ðx þ Dx; t þ DtÞ  ei ðx; t þ DtÞ þm Dt Dx eeq i ðx; tÞ  ei ðx; tÞ ¼

s

Dt 

s

eeq i ðx; tÞ  ei ðx; tÞ



ð6Þ

ð7Þ

There are different types of lattice structures employed in previous studies [18–20]. The lattice structure is represented usually by DxQy where ‘x’ denotes the number of spatial dimensions and ‘y’ indicates the number of lattice directions. The odd index after ‘Q’ in the lattice specified indicates consideration of center points of the symmetrical lattice while even index specifies lattice without center points. For example, in one dimensional lattice there can be lattices such as D1Q3, D1Q5, etc. Similarly for two dimensional problems different lattice configurations such as D2Q5, D2Q9 (Fig. 1), etc. are employed. It is observed from the inherent lattice configurations of 2-D lattices that for D2Q5 the phonons travel at a uniform speed towards neighboring nodes whereas for D2Q9 the speeds are non-uniform in order to move to the consecutive lattice nodes at the same time. Thereby weights are assigned to each lattice direction in the D2Q9 lattice to account for the non-uniformity of the lattice speeds. There are various studies available, where different weight distributions are allocated to the D2Q9 lattice [18,21–24]. In this present study a novel weight distribution is suggested. The phonons propagate simultaneously along horizontal and vertical directions along with diagonal movement for D2Q9 lattice as shown in Fig. 1. In order to calculate different weights (pi), simple mathematical analysis is employed. It can be clearly observed from geometry that the ratio of the speeds of the cross diagonal p directions in the horizontal and vertical directions is 2. Therefore taking this ratio into consideration it can be shown that p5 m5 ¼ pp6 mm63 ¼ 1, where pi and vi are the individual weights and phop m1 1

3. Numerical modeling and property estimation In this section the modeling approach and the analysis of necessary variables are explained in detail for semiconductor as well as metal. 3.1. Semiconductor material

which can be further simplified into the form as

ei ðx þ Dx; t þ DtÞ  ei ðx; tÞ ¼

Fig. 2 reflects the general structure of the LBM algorithm. Initially the material properties of the phonons such as mean free path (MFP) and relaxation time are calculated and the time step is assigned based on the grid resolution by choosing a particular lattice structure. The energy distribution function is further evaluated from a thermal boundary condition which updates itself at every time step followed by the streaming and collision steps till the convergence of the solution. The present study includes both gray model [17] and Sine dispersion model [25] and their variations are explored in terms of thermal conductivity predictions with varying porosity. The phonons can travel in any directions in three dimensional space. Under the gray model assumption (Debye model), the phonon frequency (x) is proportional to its wave number (j); expressible as

x ¼ mg j

ð8Þ

In case of silicon, the Debye model uses Eq. (8) for all frequencies and thus does not identify the difference between acoustic (LA and TA) and optical modes. Therefore the phonons have frequency independent group velocity (vg) in gray model. However in case of Sine dispersion model, the phonon dispersion relations are approximated as the vibration of a linear mono atomic chain as [26]

rffiffiffiffiffi K ja sin m 2

x¼2

ð9Þ

where K, m and a are material properties. The bulk MFP is estimated from the bulk thermal conductivity of the material and can be calculated from the kinetic theory as

3

non speeds of each lattice direction respectively. Therefore consideration of uniform phonon speeds results in the relation of the individual weights pp5 ¼ pp6 ¼ p1ffiffi2. Again it is also known that 1 3 P i pi ¼ 1 where i indicates the number of directions for a particular lattice. Here it is assumed that the center (zeroth) lattice direction does not account for phonon propagation. Thus finally it comes out that pi = 0.15 for i = 1–4 and pi = 0.1 for i = 5–8 of the D2Q9 lattice. Instead of other conventional lattices, therefore, this newly proposed weight distribution is applied to D2Q9 lattice for different two dimensional problems of various length scales.

Fig. 1. Lattice orientation: D2Q9.

481

Fig. 2. Flow chart of LBM algorithm.

482

k¼

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

1 C mK 3

ð10Þ

where C is the total heat capacity of the material at a particular temperature, m is the phonon group velocity under Debye approximation and K is the MFP, which can be represented as the product of phonon group velocity (m) and relaxation time (s). The average physical phonon and bulk properties of silicon, calculated using both Debye and Sine dispersion models at 300 K temperature are enlisted in Table 1. While solving BTE the group velocity is divided along a discrete angular resolution over the entire phase space but in case of LBM, the discretization occurs only along specified directions of the lattice. Thus in order to project the velocity vectors from a three dimensional coordinate system into the one and two dimensional planes, the estimation of effective velocity has to be exercised; bearing in mind that the calculation is exclusive for LBM. Fig. 3 represents the velocity vector in a three dimensional space, where h and / are the projections of that vector in x and y directions respectively. It can be shown that the projected velocity in y–z plane in any direction of angle / will be,

Z 1 p

ðmeff Þ2D ¼

p

m sin hdh ¼

0

2m

ð11Þ

p

Similarly for one dimension and the effective velocity will come half the actual velocity as shown below

1 2p

ðmeff Þ1D ¼

Z p

ðm sin h sin /Þ

0

Z p

sin hdhd/ ¼

0

m 2

ð12Þ

Therefore the necessary corrections should be incorporated in the velocity while calculating the MFP (K = ms) and the Knudsen number (Kn) too. The lattice equilibrium distribution in each lattice direction is related to the physical properties by the relation i

ðeeq ðx; tÞÞlat ¼ CTðx; tÞ  pi

ð13Þ

where C is the total heat capacity, T represents the physical temperature and p indicates the weight of each lattice direction. Here ‘i’ represents the index of each lattice direction ranging from 1 to m, where ‘m’ is considered as the total number of lattice directions. In case of lattices without center point ‘m’ should be replaced by ‘m  1’. Now the physical lattice temperature is calculated by

Tðx; tÞ ¼

m X

ei ðx; tÞ=C

ð14Þ

i¼1

where ei is the non-equilibrium distribution of energy distribution and C is the total heat capacity. The heat flux is also estimated by

  qn ðx; tÞ ¼ ein ðx; tÞ  eout ðx; tÞ meff cos h

ð15Þ

where superscript ‘in’ and ‘out’ describes the incoming and the outgoing energy distribution function of phonon, m is the Debye velocity, h is the angle between the phonon propagation direction and the any arbitrary direction along which the heat flux vector is going to be calculated. The average heat flux (qavg) is calculated by Table 1 Average physical bulk properties – silicon.

Phonon relaxation time (s) Phonon mean free path (K) Phonon group velocity (m) Total specific heat (C) Bulk thermal conductivity (k)

Debye model [16]

Sine dispersion model [33]

6.53  1012 s 41.8  109 m 6400 m/s 1.66  106 J/m3 K 148 W/m K

144.3  1012 s 260.4  109 m 1804 m/s 0.93  106 J/m3 K 145.6 W/m K

Fig. 3. Velocity vector in a three dimensional space.

integrating finite flux values corresponding to all spatial nodes at that plane as follows

qav g ¼

Z

nx

1

qn ðxÞdx

ð16Þ

where ‘nx’ indicates the total number of nodal points. The thermal conductivity (k) for 2-D problems can be computed by dividing the average heat flux at any plane by the temperature gradient perpendicular to the plane direction

k¼

qav g Th  Tc

ð17Þ

where ‘h’ and ‘c’ represent the leftmost hot and the rightmost cold boundaries. 3.2. Metallic material The evaluation of energy transport through porous metals requires the knowledge of property calculations of electron and phonons. Gold has been identified as the reference metal for the present study. A general expression to calculate the electron thermal conductivity (kel) of the gold is taken into consideration, which is valid over a broad range of electron temperatures [27]

 2 5=4  2  n þ 0:16 nel þ 0:44 nel kel ¼ v  el    1=2 n2el þ 0:092 n2el þ ghph

ð18Þ

k T

where, nel ¼ kBET el , hph ¼ BE ph , v = 353 W/m K and g = 0.16 for gold. In f f the lower temperature limit, ne  1 and therefore kel from Eq. (18) el reduces to the kel ¼ 315T . The heat capacity of electron is calculated T ph

from the linear relation Cel = 70Tel. The constant electron relaxation time is given by sel = 0.04 ps [28]. The values of total thermal conductivity of gold ranging from room temperature to the melting temperature are obtained from [29,30]. The temperature dependant phonon thermal conductivity is calculated to be 1% of the total thermal conductivity and the phonon specific heat is calculated based on the Debye model. The relaxation time for phonons is given by sph = 0.8 ps [31]. The electron velocity is calculated from a similar kinetic relation expressed for phonons. The frequency independent LBM for coupled electron and phonon energy distribution can be expressed as an extension of Eq. (7) mentioned below:

g i ðx þ Dx; t þ DtÞ  g i ðx; tÞ ¼ ei ðx þ Dx; t þ DtÞ  ei ðx; tÞ ¼

Dt 

sel

Dt 

sph



ð19Þ



ð20Þ

g eq i ðx; tÞ  g i ðx; tÞ eeq i ðx; tÞ  ei ðx; tÞ

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

where ‘e’ and ‘g’ symbolize phonon and electron distribution functions. The physical values of temperature, heat flux and thermal conductivity of electron can be evaluated using similar relations represented by Eqs. (14)–(17).

T ph ðx; tÞ ¼ T el ðx; tÞ ¼

m X ei ðx; tÞ=C ph i¼1 m X

g i ðx; tÞ=C el

ð21Þ ð22Þ

i¼1

qn;ph ðx; tÞ ¼ ðein ðx; tÞ  eout ðx; tÞÞmeff ;ph cos h

ð23Þ

qn;el ðx; tÞ ¼ ðg in ðx; tÞ  g out ðx; tÞÞmeff ;el cos h Z nx qav g;el ¼ qn;el ðxÞdx Z1 nx qn;ph ðxÞdx qav g;ph ¼

ð24Þ ð25Þ ð26Þ

1

qav g;el kel ¼ T h;el  T c;el qav g;ph kph ¼ T h;ph  T c;ph

ð27Þ ð28Þ

The subscripts ‘el’ and ‘ph’ indicate the physical quantities corresponding to electrons and phonons respectively. 4. Problem definition and boundary conditions The computational domain of present study consists of a regular array of uniform sized pores within the solid structure as shown in Fig. 4. However, in reality, the porous material comprises of pores of different dimensions. Also the pores may be of irregular geometries within the material. Apart from the above features, the porous material also inherently comprises of random distribution of pore networks. In order to reduce the computational intensity in terms of geometrical configuration and implementing complex boundary conditions upon the whole domain, a unit cell approach is taken into consideration. The entire domain is divided into a periodic array of single computational cells, where each cell consists of only one pore as shown in Fig. 4. In order to check the effect of different pore geometries, apart from circular pore, the results of square and rectangular pores are also compared. The porosity is calculated as the ratio of the area of pore dimension to that of the area of the square unit cell. Now in case of circular pore geometry the porosity p 2  d (/) can be calculated as / ¼ 4a2 where ‘d’ and ‘a’ the diameter the length circular pore and length of the square unit cell respectively. When ‘d’ is taken equal to ‘a’ then it can be shown that the maximum value of porosity becomes close to 0.79. However similar restrictions may not arise if we consider square and rectangular pores. Therefore, for all numerical results the porosity bandwidth is maintained up to 0.7. Knudsen number plays an important role as far as the scaling effect is concerned. Although in the present study involves othe length scales such as pore characteristic dimensions (diameter

Fig. 4. Computational domain: unit cell model with boundary conditions.

483

for circular pore – d), it is to be considered, that Kn is defined as the ratio of MFP (K) of the energy carrier to device dimension (L) unless otherwise stated. Phonon transport is examined in the unit cell using the present LBM solver employing proper boundary conditions. The heat propagation occurs from left to right; thereby periodic boundary conditions are applied along the heat flow direction and north and south boundaries are maintained as adiabatic (specular) as shown in Fig. 4. The phonon transport at the pore boundary is also implemented with suitable interfacial treatments. The interfaces are modeled as entirely diffuse or specular or a combination of both of them. Ziman proposed the following expression for calculating the specularity parameter [32]

p ¼ exp 

16p3 d2

! ð29Þ

k2

where d is the interfacial roughness and k is the phonon wavelength. The zero value of p represents a totally diffuse interface and the one value indicates the totally specular interface. A real interface is neither truly diffuse nor entirely specular, hence, the exact value of the specularity parameter of the pore boundary cannot be determined. Thus the effect of specularity parameter for a given Knudsen number has been revealed in the present study. The insulated boundary condition is maintained by applying ‘bounce forward’ scheme also known as ‘specular’ boundary condition. In this case the outgoing phonon energy distribution of a particular lattice direction equals to that of the incoming energy distribution. As far as the ‘diffuse’ boundary condition on a surface is concerned, the unknown energy distribution is calculated as an average of known energy distributions having opposite directional mobility. For example in D2Q9 lattice as depicted in Fig. 1, the specularly reflecting boundary condition is maintained at the east boundary as

e1;5;8 ¼ e2;6;7

ð30Þ

where the subscripts indicate the corresponding lattice directions. For instance while applying the diffuse boundary condition at the left east boundary the incoming energy distributions are calculated as

e1 ¼ e5 ¼ e8 ¼

e þ e þ e  2 6 7 3

ð31Þ

While implementing periodic boundary condition at the west boundary

  eq e1;5;8 ðxÞ ¼ e1;5;8 ðnxÞ þ eeq 1;5;8 ðxÞ  e1;5;8 ðnxÞ

ð32Þ

where x and nx denotes the initial and final spatial nodes and 1, 5 and 8 indicates the lattice directions. The energy distribution for neither completely specular nor diffuse surface can be estimated using the following expression:

e1 ¼ p  e2 þ ð1  pÞ 

e þ e þ e  2 6 7 3

ð33Þ

where the subscripts indicate the lattice directions and p represents the specularity parameter. The pore interface can be regarded as purely specular, diffuse or neither of them based upon the values of the specularity parameter (p). The maximum p value of 1 indicates the specular surface indicating very smooth surface with negligible roughness (low d values in Eq. (19)), while a zero value of p stands for to the diffuse one thereby representing very rough surface (high d values in Eq. (19)). Any intermediate value between zero and one of p correspond the interface being neither entirely diffuse nor specular, which is observed in practice (intermediate d values correspondiing to moderate surface roughness in Eq. (19)). It can be clearly observed that assigning the p value as one or zero in the Eq. (23) leads to either of the Eqs. (30) and (31).

484

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

5. Results and discussion In this section the present LBM model is validated against other computational and experimental data available in the literature. The entire study has been conducted using the present LBM solver from a Kn regime of 0.01–0.4 [33]. The role of different pore configurations and their performances are compared for a range of Knudsen numbers (Kn = 0.1–0.4) using both Debye and Sine dispersion models. Besides to understand the electron–phonon interaction, especially at sub continuum, energy transport of porous metallic films is also exercised. Silicon is taken as the reference material for semiconductor micro-nano porous structure, where the phonon effects are dominant, while gold is considered for the metallic study. 5.1. Semiconducting nano porous structure – silicon To validate our methodology a 2-D steady state phonon transport is solved within the unit cell domain. The phonon transport is solved using the conventional D2Q9 lattice combined with newly proposed weight configuration. The heat is transported from left to right plane as shown in Fig. 4 maintained at initial temperatures of 301 K and 300 K respectively. Initially the grid independence study has been performed with four spatial choices, starting from 100  100, 400  400, 800  800 and 1000  1000 based on estimated thermal conductivity. The corrspending thermal conductivities observed are 162 W/m K, 155 W/m K, 152 W/m K and 148.2 W/m K. As the bulk thermal conductivity of silicon is close to 148 W/m K value, 1000  1000 is selected as the optimum grid resolution. Also it was determined that 1000  1000 grid points would adequately resolve the temperature profile. 5.1.1. Validation Initially a square pore geometry is considered within the square unit cell to compare with the previous data. The LBM is compared with the results computed by the Chung et al. [7] for pore size of 1 nm and 3 nm using BTE. It can be clearly observed from the Fig. 5 that Chung et al. [7] had solved the phonon transport using fundamental BTE incorporating the pore network randomness. Gesele et al. [3] proposed the following expression for calculating the effective thermal conductivity to integrate the pore network randomness

keff ¼ kð1  /Þ3

ð34Þ

where / is the porosity and k is the thermal conductivity calculated from numerical model. The reason for integrating this expression with the numerical predicted values is that the inaccuracy in predicting the thermal conductivity can be minimized while calculating the thermal properties of submicron porous structures involving non homogeneous random pore distribution. Besides, it can be also mentioned that the computed BTE results show a marked improvement while incorporating this relation as observed by Chung et al. [7]. The LBM results reflected in Fig. 5 also include the pore network randomness proposed by Gesele et al. [3]. However for low values of porosity this expression may not able to predict the thermal conductivity accurately. The predicted thermal conductivity plots are compared with the existing experimental values available in the literature. The detailed parameters corresponding to each of these experimental observations indicated in Fig. 5 are enlisted in Table 2. The experimental values are considered on the basis of pore dimensions where the pore sizes are close to around 10 nm and below thereafter. When the predicted values from the present LBM solver are compared with the measured values of Benedetto et al. [2] and Drost et al. [4], a general agreement is observed among them. How-

Fig. 5. Comparison of LBM with BTE [6] and experimental observations [1–3].

ever Fig. 5 also reflects that the estimated thermal conductivity values of the LBM model over predicts by one order of magnitude than the experimental observations of Gesele et al. [3], especially for high porosity values. In fact LBM exhibits a decent agreement with the BTE values considering pore random network [7] except for higher porosity values greater than 0.7. Therefore it can be concluded that LBM can successfully predict the general trend of size effects in terms of porosity variations. 5.1.2. Effect of pore interface While comparing the size effects for various choices of specularity parameters (p = 0.0, 0.5, 1.0), it can be observed that from Fig. 6 that for a lower Knudsen number of 0.1 distinct profiles of thermal conductivity variations are noticed. The maximum value of thermal conductivities is reported while the pore interface is considered as totally specular (p = 1.0) and for any intermediate choice of p = 0.5, the prediction remains in between these extreme limits. However as the pore dimensions approach towards subcontinuum it can be pointed out that the effect of the specularity parameter has negligible influence in estimating the thermal conductivity values. The collapse of all the three p values (0.0, 0.5, 1.0) into a single profile in case of Kn of 0.4 does yield the impression of the strong ballistic effects, which may not able to distinguish regarding the scattering information provided at the pore interface, particularly at high Knudsen numbers. Fig. 6 also reveals that with increasing porosity values the thermal conductivity reduces consistently. 5.1.3. Comaprative study of various pore configurations Fig. 7 reflects the comparative study of both square and circular pore geometries for two different Knudsen numbers of 0.1 and 0.4 with varying specularity parameters. It can be noticed that for Kn of 0.1 the predictions of the thermal conductivity in case of circular

Table 2 Experimental data for Fig. 5. Samples

A

B

I

a

b

c

d

Porosity (/) Pore size (nm) k (W/m K)

0.6 >12

0.5 12

0.4 3

0.64 9.0 ± 3.0

0.64 1.7 ± 0.5

0.71 2.0 ± 0.3

0.79 2.7 ± 0.3

2.5

3.9

1.2

0.79

0.196

0.135

Experimental observations: [3].

A, B

Bendetto et al. [2], IDrost et al. [4],

0.058 a,b,c,d

Gesele et al.

485

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

Fig. 6. Comparison different specularity parameters (0.0, 0.5and 1.0) for various Knudsen numbers (0.1, 0.2 and 0.4).

pore also reflect similar discrete patterns with that of the square one. However the numerical values of the circular ones anticipate higher values than that of its compatriot and the deviation between them increases with increasing values of porosity. Besides for higher Kn of 0.4, similar observations are noticed for the circular ones compared to the square pore, however with higher values and likewise divergence between the two for high porosity. In order to understand the reason behind this it has been attempted to compare both the circular and square pore configuration with the rectangular one. The case of Kn of 0.1 having porosity of 0.4 with entirely specular pore interface boundary is chosen for the comparison among pore configurations. The predicted LBM results are enlisted in Table 3. The comparative study unveils that for a given porosity the characteristic length scale of the pore is higher in case of circular one than the square one. Further, while employing the rectangular pore dimensions it can be noticed by varying the length of the sides a lower value of thermal conductivity is achieved maintaining the same conditions. In the Table 3, for the rectangular configuration ‘a’ and ‘b’ notations indicate the base and height of a rectangular

pore dimension. It is also observed that when the height of the rectangular pore is changed from 140 nm to 226 nm maintaining the same porosity for a particular condition, thermal conductivity value of the latter becomes half of the former. Therefore in order to achieve a porous material of low thermal conductivity such customization of pore geometry can expose a novel direction in tailoring material properties accordingly. In order to understand the reason of the reduction of thermal conductivity from circular geometry to the rectangular one, the physical phenomena of phonon transport happening within the computational domain needs interpretation as well as an explanation. The phonons travel along the temperature gradient from the left boundary towards the right as shown in Fig. 4. The phonons scatter with the interfacial pore boundary on their way towards the temperature gradient. Actually the phonons can travel in any directions. However while solving the phonon transport in LBM, phonons are guided only along the adjacent lattice directions, which for the D2Q9 lattice are either horizontal, vertical or crossdiagonal directions. For further analysis only the horizontal movements of phonons are considered along the temperature gradient directions. Similar analogy can be drawn for other directional movements. It can be observed from the Fig. 8 that when the phonons are excited within the implemented temperature gradient across the boundary, the pore interface region acts as an impediment in the process of energy transport. Therefore considering only horizontal phonon propagation as mentioned in the above section, phonons apart from those emitting out of the projected plane of ‘ae’, having east bound movement will not scatter at the pore boundary, rather they will exchange energy with the right end boundary of the computational domain. This phenomena hold good for all the three geometric configurations (circle, square and rectangle) considered here. Further it is already observed from Table 3 that for rectangular pore, customizing the height of the pore can lead to low thermal conductivity values for the same computational domain. In fact it can be observed that the rectangular porous structures yield lower thermal conductivity values than that of the circular ones of the same porosity. It is interesting to note that for the prescribed two dimensional transport along the temperature gradient, the effect of scattering area becomes a function of interfacial pore length scale such as diameter for circular one, length for square and height for rectangular pore. It is noteworthy that for a few selected cases of the rectangular pore, the height of the pore is more than that of the diameter of the circular pore; whereas in other cases the height of the rectangular pore is less than the circular pore diameter. It is explained what factors influence the lowering of the thermal conductivity when the height of the rectangular pore is variable with respect to the diameter of the circular pore. It can be reasoned that for greater height the rectangular one, more number of phonons encounter the rectangular pore interface

Table 3 LBM results for Kn = 0.1, / = 0.4 and p = 0.0. Circular

D = 202 nm k = 120.1 W/m K

Fig. 7. Comparison of square and circular pore configurations for different Knudsen numbers of 0.1 and 0.4 with varying specularity parameters (0.0, 0.5 and 1.0).

Square

A = 180 nm k = 78.4 W/m K

Rectangular A  B (nm  nm)

k (W/m K)

140  226 148  212 154  204 170  188 188  170 204  154 212  148 226  140

52.3 60.8 64.8 74 84.2 92.2 96.3 100.9

Note: D – diameter of a circle; A – length of a side of square; A, B – lengths of the sides of a rectangle; k – thermal conductivity.

486

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

Fig. 8. Different pore configurations: (i) circle, (ii) square and (iii) rectangle.

and scatter to and fro between the pore interface and the left domain boundary. Therefore due to increase in the length of the projected scattering plane ‘ae’, the rectangular one will experience more number of phonon scatters at the pore boundary than the circular one. Hence the increasing scattering area on the pore interface plane contributes in reducing the thermal conductivity of the porous material as shown in Fig. 9. Therefore, it is observed that the available scattering area of the pore becomes a key criteria in scaling down the thermal conductivities and the predicted thermal conductivity varies inversely with the scattering area for rectangular and square geometry. Further questions can be raised, how the thermal conductivity shows lower value when the height of the rectangular pore is less than the circular pore diameter. For instance in case of circular porous geometry, the phonons propagating along ‘cd’ plane will suffer more number of scattering events than the ‘ab’ plane, even though both of them are confronted by the same pore interface. Therefore for circular pore geometry the phonons emitted from the ‘ae’ plane will experience non uniform scattering in terms of its frequency. However cases like square and rectangular pores do not feature similar phenomena due to their inherent geometrical configuration as the pore interface plane is parallel to the plane of phonon emission. For both of them the scattering rate is uniform and the frequency of the scattering event remains the same unlike the circular one. Therefore, for a constrained porosity; the square and rectangular pores experience more frequent scattering phenomena which ultimately reduces their thermal conductivity than the circular geometry.

5.1.4. Comparison between Debye and Sine dispersion models In addition to the Debye model explained earlier, it has been attempted to compare the same with the Sine dispersion model to check how each of their behavior influence the solutions while predicting thermal conductivity. It can be observed from the Fig. 10 that the Sine dispersion model have insignificant variations from the Debye model when the thermal conductivities of two different Knudsen numbers (0.1 and 0.4) are compared with the varying specularity parameter. Although for low Kn of 0.1 the profiles have distinct nature, there remains no difference when the thermal conductivities are plotted for Kn of 0.4 against porosity. Hence it can be concluded that both of the two models function in a similar manner while calculating the thermal properties. Thus, when the length scale of the system approaches sub micron limits, the role of specularity parameter on pore interface scattering becomes insignificant due to the ballistic nature of phonon transport. However it is to be mentioned that when the Sine dispersion model is implemented, the selection of pore geometry still remains a central criteria in controlling the thermal conductivity as reflected in Fig. 7.

Fig. 9. Variation of thermal conductivity for different scattering areas of the rectangular pore for Kn of 0.1 and porosity of 0.4.

Fig. 10. Comparison of Debye and Sine Dispersion models for Knudsen numbers of 0.1 and 0.4 for three different specularity parameters (0.0, 0.5 and 1.0).

5.2. Metallic nanostructures – gold The energy transport within metallic films of submicron scale is different from that of the semiconductor one. For semiconductor films, the phonon play a dominant role in energy transfer, whereas for metals electron behave as the principal energy carrier. The successful performance of the present LBM solver while solving

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

phonon transport in the continuum is already established and demonstrated by Chattopadhyay et al. [33]. Fig. 11 shows the non dimensional temperature profile of electron transport in bulk gold solid structure. The linear trend of this plot reflects the continuum behavior as expected. In order to calculate the thermal properties of gold, the predicted value of thermal conductivity of electron part comes out as 310 W/m K, which is very close to the bulk thermal conductivity of gold 318 W/m K. The electronic contribution of thermal conductivity of gold is almost 100 times higher than the phonon contribution at room temperature for bulk case [34]. Considering this fact it can be inferred that the present LBM solver can capture the electron transport and the same model can be extended to study the combined contribution of both electrons and phonons towards energy transport. It has been observed that electrons have a major share in the thermal conductivity estimation at the bulk regime. However when the device dimension approach the sub continuum limit the complementary role of phonons cannot be neglected. It would be therefore significant to understand the scaling of the properties of each of them in the energy transport process. The electron and phonon velocity within the gold can be calculated using the property calculations mentioned earlier. It is to be observed that electrons have almost three orders of magnitude higher velocity than its counterpart (vel  1  106 m/s and vph  2.2  103 m/s). The mean free paths of electron and phonon also reflect considerable differences (Kel  41.6 nm and Kph  1.75 nm). Knudsen number in this section is calculated based on the mean free path of the electron as a reference for further analysis. Hence it can be understood that for a given Knudsen number of 0.1, the reference length scale assumes closes to 400 nm. Thus in case of Kn of 0.1 the electron may reflect size effects as it is already advancing towards ballistic limits, however, at the same length scale the phonon is well within the continuum regime. The numerical examples are solved same as aforementioned with different pore configurations for various Knudsen numbers from 0.1 to 0.4 by changing the specularity parameter. Fig. 12 shows a reference case of porosity 0.4 having square pore, where the thermal conductivity contribution of each energy carrier of three different Knudsen numbers (Kn = 0.1, 0.2, 0.4) is represented in the comparison bar chart. It can be noticed that the phonon share in thermal conductivity prediction remains unaltered for all the Knudsen numbers. However the electron parts indicate a considerable decrease in the thermal conductivity value. In fact

Fig. 11. The bulk temperature profile of electron transport in gold.

487

Fig. 12. Comparison of electron and phonon thermal conductivities for various Knudsen numbers (0.1, 0.2 and 0.4) in the case of square pore having porosity of 0.4.

the estimated thermal conductivity from electronic component decrease to approximately to the half when the Knudsen number is varied from 0.1 to 0.4. The reason for the gradual reduction of the electron thermal conductivity is due to size effects. As the characteristic length starts diminishing, the ballistic nature of electron transport takes place, thereby reducing the thermal conductivity value. However even if the Kn changes from 0.1 (400 nm) to 0.4 (100 nm), phonons recognize the length scales within the continuum limit as the phonon MFP is more than twenty times lower than the electron, as discussed above. Hence it can be concluded that reduction of Kn from 0.1 to 0.4 has negligible impact upon calculating the phonon thermal conductivity values, whereas the electron suffers a considerable reduction thereby even decreasing the total thermal conductivity of porous gold film. In order to study the effects of various pore configurations upon the electron–phonon combined energy transport, the square and circular pore geometries are taken as reference where the porosity is considered as 0.7. It is observed from Fig. 13 that for Kn of 0.1, the predicted electron thermal conductivity reduces to an approximate 2.5 times from 249 W/m K to 100 W/m K, when the pore geometry is altered from circular to the square one. Also for Kn of 0.4, the decrease in predicting the thermal conductivity of the metal compared to Kn of 0.1 exhibits the size effects clearly which is consistent irrespective the choice of pore geometry. Further when the phonon thermal conductivities of the circular and the square pores are compared, it is noticed that there is a minor reduction in the predicted values which ensures the geometric effect. Therefore it is confirmed that the suitable pore configurations influence the solutions in estimating the thermal conductivity of the porous materials. Moreover it is noteworthy to detect that for low values of Knudsen number (Kn of 0.1), the ratio of the electron to phonon thermal conductivity value is quite high, thereby conforming the dominant electron transport compared to that of the phonon. However with the inception of size effcts

Fig. 13. Comparison of electon and phonon thermal conductivities for Knudsen numbers of 0.1 and 0.4 in the case of square and circular pore having porosity of 0.7.

488

A. Chattopadhyay, A. Pattamatta / International Journal of Heat and Mass Transfer 72 (2014) 479–488

at higher Knudsen numbers, the predicted thermal conductivities of them become comparable by monitoring the pore configuration. Hence forth it can be concluded that as the device dimension approach ballistic limits and for large values of porosity, the phonon contribution towards thermal conductivity cannot be neglected compared to the electron contribution. 6. Summary and conclusions The study focuses upon a detailed comparison of different nanopore configurations with varying interfacial pore scattering value for modeling two dimensional heat transfer using Lattice Boltzmann Mthod (LBM). It is found that LBM can predict the thermal conductivity of semiconductor porous structures such as silicon with reasonable accuracy at sub continuum regime when compared to experimental values. The estimated thermal conductivities, calculated with various values of the specularity parameter at the pore boundary for a particular geometry show a distinct nature for low Knudsen number cases. However when the length scale of the system approaches sub micron limits, the role of specularity parameter on pore interface scattering becomes insignificant due to the ballistic nature of phonon transport. When the squared pore configuration in the unit cell model is replaced by circular pore with same porosity, it is observed that the predicted thermal conductivity of circular one results in reasonably higher value than that of its counterpart. Besides while the pore geometry is transformed into a rectangle it is also noticed that even lower thermal conductivity than the square pore can be achieved due to an increase in the number of phonon scattering events with higher frequency as the scattering area varies accordingly. It is also noteworthy that for the same porosity the estimated thermal conductivity varies inversely with the scattering area. As far as the variation of model is concerned, it is observed that for low Knudsen numbers the predicted thermal conductivity profiles of the Sine dispersion model indicate a marginal departure from that of the Debye model. Further when the energy transport is solved in porous metallic structure, the electron transport shows considerable size effects whereas phonons have minimal interference in terms of estimated thermal conductivity, due to their variance in individual MFPs. Thereby this study emphasizes that, the contribution of phonon towards thermal conductivity in metals cannot be neglected at sub-continuum having large values of porosity. Finally it can be infered that porosity of structure play a dominant role in scaling down the thermal conductivity. However for a fixed porosity the variation of the pore configurations of a micro-nano structure aids in tailor making the thermal properties, irrespective of the choices of materials (for both metal and semiconductor); which promises a new way of energy conversion mechanism. Acknowledgements The authors would like to acknowledge the resources of the VIRGO supercomputing cluster at the High Performance Computing Centre (HPCE), IIT Madras for the simulations. References [1] R. Venkatasubramanian, Recent Trend in Thermoelectric Materials Research III: Semiconductors and Semimetals, Academic Press, San Diego, 2001.

[2] G. Bendetto, L. Boarino, R. Spagnolo, Evaluation of thermal conductivity of porous silicon layers by a photo acoustic methods, Appl. Phys. A 64 (1996) 155–159. [3] G. Gesele, J. Linsmeir, V. Drach, J. Fricke, R. Arens-Fischer, Temperature dependent thermal conductivity of porous Silicon, J. Phys. D: Appl. Phys. 30 (1997) 2911–2916. [4] A. Drost, P. Steiner, H. Moser, W. Lang, Thermal conductivity of porous silicon, Sens. Mater. 7 (1995) 111–120. [5] A. Wolf, R. Brendel, Thermal conductivity of sintered porous silicon films, Thin Solid Films 513 (2006) 385–390. [6] D. Song, G. Chen, Thermal conductivity of periodic microporous silicon films, Appl. Phys. Lett. 84 (2004) 687. [7] J.D. Chung, M. Kaviany, Effects of phonon pore scattering and pore randomness on effective conductivity of porous silicon, Int. J. Heat Mass Transfer 43 (2000) 521–538. [8] J.-H. Lee, J.H. Grossman, J. Reed, G. Galli, Lattice thermal conductivity of nanoporous Si: molecular dynamics study, Appl. Phys. Lett. 91 (2007) 223110. [9] Q. Hao, G. Chen, M.S. Jeng, Frequency-dependent Monte Carlo simulations of phonon transport in two-dimensional porous silicon with aligned pores, J. Appl. Phys. 106 (2009) 114321. [10] J. Randrianalisoa, D. Ballis, Monte Carlo simulation of cross-plane thermal conductivity of nanostructured porous silicon films, J. Appl. Phys. 103 (2008) 053502. [11] R. Prasher, Thermal conductivity of composites of aligned nanoscale and microscale wires and pores, J. Appl. Phys. 100 (2006) 034307. [12] R. Prasher, Transverse thermal conductivity of porous materials made from aligned nano- and microcylindrical pores, J. Appl. Phys. 100 (2006) 064302. [13] S.S. Ghai, W.T. Kim, R.A. Escobar, C.H. Amon, A novel heat transfer model and its application to information storage systems, J. Appl. Phys. 97 (2004). 10P703 (1–3). [14] R.A. Escobar, S.S. Ghai, M.S. Jhon, C.H. Amon, Multi-length and time scale thermal transport using the Lattice Boltzmann method with application to electronics cooling, Int. J. Heat Mass Transfer 49 (2006) 97–107. [15] P. Heino, Lattice-Boltzmann finite difference model with optical phonons for nanoscale thermal conduction, Comput. Math. Appl. 59 (2010) 2351–2359. [16] A. Christensen, S. Graham, Multiscale Lattice Boltzmann modeling of phonon transport in crystalline semiconductor materials, Numer. Heat Transfer B: Fundam. 57 (2010) 89–109. [17] A. Nabovati, D.P. Sellan, C.H. Amon, On the Lattice Boltzmann method for phonon transport, J. Comput. Phys. 230 (2011) 5864–5876. [18] Y.H. Qian, D. D’Humieres, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (6) (1992) 479–484. [19] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [20] J. Zhang, Lattice Boltzmann method for microfluidics: models and applications, Microfluid Nanofluid 10 (2011) 1–28. [21] A.A. Mohammad, Lattice Boltzmann Method, Springer, London, 2011. [22] P.L. Bhatnagar, E.P. Gross, M. Krook, A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems, Phys. Rev. 94 (1954) 511–525. [23] M.E. Kutay, A.H. Aydilek, E. Masad, Laboratory validation of Lattice Boltzmann method for modeling pore-scale flow in granular materials, Comput. Geotech. 33 (2006) 381–395. [24] S.S. Chikatamarla, I.V. Karlin, Lattices for the Lattice Boltzmann method, Phys. Rev. E 79 (1–18) (2009) 046701. [25] G. Chen, Size and interface effects on thermal conductivty of superlattices and periodic thin-film structures, J. Heat Transfer 119 (1997) 220–229. [26] G. Chen, Nanoscale Energy Transport and Conversion, Oxford University Press, New York, 2005. [27] S.I. Anisimov, B. Rethfeld, On the theory of ultrashort laser pulse interaction with metal, SPIE 3093 (1997) 192–203. [28] T.Q. Qiu, C.L. Tien, Heat transfer mechanisms during short-pulse laser heating of metals, ASME J. Heat Transfer 115 (1993) 835–841. [29] Y.S. Touloukian, R.W. Powell, C.Y. Ho, P.G. Klemens, Thermal conductivity, Thermophys. Prop. Matter 1 (1970) 132–137. [30] Y.S. Touloukian, R.W. Powell, C.Y. Ho, P.G. Klemens, Thermal conductivity, Thermophys. Prop. Matter 1 (1970) 132–137. [31] J.K. Chen, J.E. Beraun, Numerical study of ultrashort laser pulse interactions with metal films, Numer. Heat Transfer A 40 (2001) 1–20. [32] J. Ziman, Electrons and Phonons, Oxford University Press, London, 1960. [33] A. Chattopadhyay, A. Pattamatta, Estimation of an appropriate lattice structure for phonon transport using Lattice Boltzmann method, in: Proceedings of ASME Summer Heat Transfer Conference, Paper HT2013-17188, Minneapolis, MN, USA, 2013. [34] A. Pattamatta, C. Madnia, A comparative study of two-temperature and Boltzmann transport models for electron–phonon nonequilibrium, Numer. Heat Transfer: A 55 (2009) 611–633.