Study on Lattice Thermal Conductivity of Silicon Thin Film with Aligned Nano-pores

Available online at www.sciencedirect.com

ScienceDirect Energy Procedia 105 (2017) 4915 – 4920

The 8th International Conference on Applied Energy – ICAE2016

Study on lattice thermal conductivity of silicon thin film with aligned nano-pores Qi Liang, Ya-Ling He*, Yi-Peng Zhou, Tao Xie Key Laboratory of Thermo-Fluid Science and Engineering of MOE, School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China

Abstract A frequency-dependent phonon Boltzmann transport equation solver is presented to numerically study the lattice thermal conductivity of nano-porous silicon thin film based on the Discrete Ordinates Method. We find that not only characteristic size and porosity, but also the shapes of both unit cell and pore do affect the lattice thermal conductivity. Therefore, we arranged a series of computational cases by the orthogonal design method to investigate the influence of geometric parameters Lx, Ly, ax/Lx and ay/Ly (Lx is unit cell length, Ly unit cell width, ax pore length and ay pore width). Furthermore, a non-linear regression model is established depending on the data obtained from those computational cases. The result shows that lattice thermal conductivity of nano-porous silicon thin film decreases obviously with the increase of ax/Lx and ay/Ly. Among these four geometric parameters, ay/Ly is the most significant factor while Lx and Ly have little effect on lattice thermal conductivity. The proposed regression model can offer useful suggestion for the fabrication of nano-porous silicon thin film with lower lattice thermal conductivity. ©©2017 Published by Elsevier Ltd. This an open access 2016 The Authors. Published byisElsevier Ltd. article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Selection and/or peer-review under responsibility of ICAE Peer-review under responsibility of the scientific committee of the 8th International Conference on Applied Energy.

Keywords: Lattice thermal conductivity, discrete ordinates method, nano-porous thin film, orthogonal design method.

1. Introduction Suffering from the increasingly energy crisis and the resulting environment pollution, the usage of renewable energy, especially efficient utilization of solar energy such as photovoltaic cells, has aroused great interest. To further improve the efficiency of photovoltaic cells in harvesting and converting solar energy, several hybrid systems combining photovoltaic (PV) cell and thermoelectric (TE) module were promoted recently [1][2]. Researches showed that better thermoelectric materials would allow for an efficiency increase of up to 50%[3]. The ideal efficiency of the TE module is determined by the

* Corresponding author. Tel.: +86-29-8266-5930; fax: +86-29-8266-5445. E-mail address: [email protected].

1876-6102 © 2017 Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Peer-review under responsibility of the scientific committee of the 8th International Conference on Applied Energy. doi:10.1016/j.egypro.2017.03.981

4916

Qi Liang et al. / Energy Procedia 105 (2017) 4915 – 4920

dimensionless figure of merit defined as ZT=(σS2/λ)T, where σ is the electrical conductivity, S the Seebeck coefficient, λ the thermal conductivity, and T is the temperature[4]. Si is the most commonly used semiconductor material, which is abundant and well-engineered, but has little ZT for its high thermal conductivity. Recent work has shown that thin-film nanoengineering may significantly decrease the thermal conductivity of silicon with slight effect on the electrical conductivity and Seebeck coefficient due to different size effects between phonons and electrons[5][6][7]. To investigate the phonon transport characteristics in the nanostructured materials and further decrease lattice thermal conductivity, various simulation methods, such as lattice Boltzmann method[8], molecular dynamics[9], Monte Carlo (MC) method[10][11], Boltzmann transport equation[12][13] etc., have been applied. In this paper, the frequency-dependent phonon Boltzmann transport equation (PBTE) is solved based on the Discrete Ordinates Method (DOM) to reveal the influence of pore structure and unit cell structure on phonon transport property in the nano-porous silicon thin film. 2. Model and methodology description 2.1. Geometric model The structure of silicon thin film with aligned pores is given in Fig. 1. The direction perpendicular to the surface of thin film is called cross-plane direction, while the direction parallel to the surface is inplane direction. In this work, only in-plane heat transfer is investigated and the corresponding unit cell is shown in Fig. 1 (b). The heat is assumed to flow along the positive direction of x axis, which means that temperature at the left-side boundary is higher than that at the right-side boundary. The length L of the unit cell is the characteristic length, which denotes periodic size of the porous structure, and a is the length of pores. The porosity of the porous structure is defined as Φ=a2/L2. 2.2. Description of the numerical method By defining the phonon intensity I  p, Z, r , sˆ f Q  p, k Z  p, k D Z 4 , the PBTE is transformed into the equation of phonon radiative transport (EPRT). Here, h 2π is the reduced Planck constant, Q  p, k the group velocity of phonon with phonon branch p and wave vector k, Z  p, k the angular frequency of phonons, D Z the density of states for phonons, and as shown in Fig. 2, sˆ [ ,K , P  cosT ,sin T cos\ ,sin T sin\ is the unit vector of phonon transport direction, where θ is the polar angle and ψ is the azimuthal angle. In the present study, only steady state heat conduction is considered and the corresponding three dimensional EPRT can be written as I0 wI wI wI I [ K  P  (1) wx wy wz Q  p, k W  p, Z , T Q  p, k W  p, Z , T , where W  p, Z , T is the phonon relaxation time, I0 the local equilibrium phonon intensity when f is the equilibrium distribution following the Bose-Einstein distribution[14]. To solve this equation, the Discrete Ordinates Method (DOM) is employed here, which is a method widely used for solving the Thermal Radiative Transport Equation (TRTE)[15]. In the Discrete Ordinates Method, Eq. (1) is solved for a set of n different directions sˆi , that is wI wI wI Ii I0 [ i i  Ki i  Pi i  ˈi 1, 2, , n , (2) wx wy wz Q  p, k W  p, Z , T Q  p, k W  p, Z , T where Ii

I  p, Z, r , sˆi indicates the phonon intensity for direction ŝi. Frequency-dependent equilibrium

phonon intensity is defined as I0  p, Z, r

¦

n i 1

wi I  p, Z, r, sˆi 4π , where wi is the quadrature weight

4917

Qi Liang et al. / Energy Procedia 105 (2017) 4915 – 4920

for direction ŝi. The total equilibrium phonon intensity is obtained by integrating the frequency-dependent equilibrium phonon intensity over phonon frequency and branch. The effective thermal conductivity is calculated by O = qx L 'T , in which qx is the average heat flux along x direction. y K sˆ \

0

T

x

[

P

Fig. 1 The sketch and geometric model for nano-porous silicon thin film. (a) Schematic of

Fig. 2. Coordinate system used

silicon thin film with aligned nano-pores. (b) Geometric model for In-plane condition.

for phonon transport simulation.

In our work, S8-approxiamtion is adopted, which has been proved adequate for the simulation of phonon transport in silicon films[12]. And diffuse reflection boundary condition is employed for pore boundaries, while the boundaries of the unit cell are considered as period boundary condition[16]. For silicon, phonon dispersion along the (100) lattice direction is used here, assuming that the phonon transport is isotropic in silicon. Among the six phonon branches, we only take three acoustic phonon branches in consideration, including two transverse acoustic (TA) branches and one longitudinal acoustic (LA) branch, duo to the optical phonons contribute little to the thermal conductivity at room temperature. The frequency of phonons for silicon is calculated by Pop’s relation[17], written as Z Z0 Q s k  ck 2 , where the parameters ω0, νs and c are fitting coefficients given in Table 1. The maximum wave vector in the first Brillouin zone is kmax=2π/a, and a=5.43 Å is the lattice constant for Silicon. For internal scattering, only the impurity scattering and phonon-phonon scattering are considered and the charge carrier scattering is ignored for the assumption of light doping. The relaxation time of impurity scattering is expressed as W i1 AiZ 4 , where Ai=1.32 h 10-45 s3[18]. And relaxation times of phonon-phonon scattering for Normal process and Umklapp process are calculated by the combined relaxation time 1 W N1  W U1 . Based on molecular dynamics simulation results obtained by Henry et W NU ,.defined as W NU 2 al.[19], W NU ANU Z 2π T  b , where ANU and b are given in Table 1. Then, the overall phonon relaxation time is attained by W 1 W i1  W N1  W U1 according to the Matthiessen’s rule. Table 1. Fitting coefficients for phonon dispersion and relaxation time of Silicon

TA branch LA branch

ω0 rad·s-1

νs m·s-1

c 10-7 m2·s-1

ANU Kb/s

b

0.00 0.00

5230 9010

-2.26 -2.00

5.07h1018 5.32h1018

1.65 1.49

3. Results and discussion 3.1. Validation of the numerical model To validate the present numerical method, the in-plane thermal conductivity of porous silicon film with porosity Φ=0.25 is calculated by our code under the circumstance of T=300 K. For comparison, the same transport properties, including the phonon dispersion and relaxation time, should be adopted as that used in Hao’s MC simulation[11]. The relaxation time used here is the same as that employed in the reference paper, while the phonon dispersion is different. Hence, the phonon dispersion needs to be

4918

Qi Liang et al. / Energy Procedia 105 (2017) 4915 – 4920

modified. Hsieh et al.[16] has fitted the dispersion with quadratic polynomial in the form given above. For the TA branch, ω0, νs and c are 0 rad/s, 6.14h103 m/s and -1.64h10-7 m2/s, respectively. For the LA branch, ω0, νs and c are 0 rad/s, 8.97h103 m/s and -1.37h10-7 m2/s, respectively. The results predicted by our code is plotted in Fig. 3, as well as those obtained by MC simulation and by DOM simulation. From the diagram, we find that present predictions are in good agreement with the results of Hao’s MC simulation and Hsieh’s DOM simulation. 3.2. Lattice thermal conductivity of nano-porous silicon thin film In previous researches, the shapes of unit cell and pore are usually fixed as square. To investigate the influence of the structures of unit cell and pore on phonon transport property in the nano-porous silicon thin film, the shapes are specified as rectangles instead, as shown in Fig. 4. The length (along the x axis) of unit cell is defined as Lx, while the width (along the y axis) of unit cell is Ly. And the length and width of pore are defined as ax and ay, respectively. For comparison, unit cell area and porosity are fixed. The size of unit cell is 400 u 400 nm2 for (a)-(c), 640 u 250 nm2 for (d) and 250 u 640 nm2 for (e). While pore sizes are 200 u 200 nm2, 320 u125 nm2 and 125 u 320 nm2, respectively. The simulation results show that λc<λa<λd<λe<λb, and λa, λd and λe are close to each other among the range of 31.50-32.95 W/m·K, while λc=13.62 W/m·K and λb=43.76 W/m·K. Comparing (b) with (c), large difference in lattice thermal conductivity is observed for the diverse direction of pores. For (a), (d) and (e), shapes of both unit cell and pore are different but similar lattice thermal conductivity is obtained. The results show that both unit cell structure and pore structure do affect the phonon transport property in silicon thin film. The diversion of lattice thermal conductivity between them may be aroused from different ax/Lx and ay/Ly.

O(W/m· K)

102

101 DOM, present work MC, Hao et al.[11] DOM, Hsieh et al.[15] 100 0 10

101

102 L (nm)

103

104

Fig. 3. Thermal conductivity of nano-porous silicon thin film for different unit cell lengths.

Fig. 4. Temperature distributions in silicon thin films with different unit cells and pores

3.3. Analysis based on the orthogonal design method To uncover the influence of four geometric factors (Lx, Ly, ax/Lx and ay/Ly) on the phonon transport characteristics in silicon thin film, a set of computational cases is determined based on orthogonal array L25(56). The factors and levels used are presented in Table 2. After simulations, the summed lattice thermal conductivity for each factor is obtained and plotted in Fig. 5. It shows that ay/Ly has a great influence on phonon transport and ax/Lx has a relatively large influence, while the influence of Lx and Ly are small. The summed lattice thermal conductivity decreases with the increase of ay/Ly and ax/Lx, and increases with the increase of Lx and Ly. Furthermore, the lattice thermal conductivity is fitted with a non-

4919

Qi Liang et al. / Energy Procedia 105 (2017) 4915 – 4920

linear regression model through the non-linear least-squares method, as shown in Eq. (3). The corresponding R-square is 0.98, indicating that the obtained regression model is reliable. 0.768 0.134 O ª« 2.033L0.004  ax Lx  a y Ly º¼» 0.021  e Ly  93.845 . x (3) ¬





Table 2. Factors and levels used in this study

Code

Factors

A B C D

Lx(nm) Ly(nm) ax/Lx ay/Ly

Level 1

2

3

4

5

400 400 0.1 0.1

800 800 0.3 0.3

1200 1200 0.5 0.5

1600 1600 0.7 0.7

2000 2000 0.9 0.9

Osum (W/m· K)

400

300

200

100

0

400

800

1200 1600 Lx (nm)

2000

400

800

1200 1600 Ly (nm)

2000

0.1

0.3

0.5 ax/Lx

0.7

0.9

0.1

0.3

0.5 ay/Ly

0.7

0.9

Fig. 5. Summed lattice thermal conductivity for each factor

4. Conclusions In this paper, the frequency-dependent phonon Boltzmann transport equation is solved for study of phonon transport characteristics in silicon thin film using the Discrete Ordinates Method. Depending on the PBTE solver, both unit cell structure and pore structure are studied based on the computational cases arranged by orthogonal design method. The result shows that when the unit cell size is fixed, the shape and direction of pores have great influence on phonon transport property. The lattice thermal conductivity of silicon thin film decreases with the increase of ax/Lx and ay/Ly, and increases with the increase of Lx and Ly. Among these four geometric parameters, ay/Ly is the major factor while Lx and Ly have little effects on the phonon transport. Finally, a non-linear regression model is fitted for lattice thermal conductivity of silicon thin film with aligned nano-pores, which offers useful suggestion for the fabrication of nanoporous silicon thin film with lower lattice thermal conductivity to obtain better thermoelectric material. Acknowledgements This work is supported by Major Program of the National Natural Science Foundation of China (Grant No. 51590902). References [1] Hsueh TJ, Shieh JM, Yeh YM. Hybrid Cd-free CIGS solar cell/TEG device with ZnO nanowires. Prog Photovolt Res Appl 2015;23:507-512.

4920

Qi Liang et al. / Energy Procedia 105 (2017) 4915 – 4920

[2] Wang Y, Su S, Liu T, et al. Performance evaluation and parametric optimum design of an updated thermionic-thermoelectric generator hybrid system. Energy 2015;90:1575-1583. [3] Van Sark W. Feasibility of photovoltaic-thermoelectric hybrid modules. Appl Energ 2011;88:2785-2790. [4] Goldsmid, HJ. Introduction to thermoelectricity. Vol. 121. Springer Science & Business Media, 2009. [5] Tang J, Wang HT, Lee DH, et al. Holey silicon as an efficient thermoelectric material. Nano Lett 2010;10:4279-4283. [6] Yu JK, Mitrovic S, Tham D, et al. Reduction of thermal conductivity in phononic nanomesh structures. Nat nanotechnol 2010;5:718-721. [7] Hopkins PE, Reinke CM, Su MF, et al. Reduction in the Thermal Conductivity of Single Crystalline Silicon by Phononic Crystal Patterning. Nano Lett 2011;11:107-112. [8] Sellan DP, Turney JE, Mcgaughey AJH, et al. Cross-plane phonon transport in thin films. J Appl Phys 2010;108:113524. [9] Chalopin Y, Esfarjani K, Henry A, et al. Thermal interface conductance in Si/Ge superlattices by equilibrium molecular dynamics. Phys Rev B 2012;85:195302. [10] Jeng MS, Yang R, Song D, et al. Modeling the Thermal Conductivity and Phonon Transport in Nanoparticle Composites Using Monte Carlo Simulation. J Heat Trans-T ASME 2008;130:042410. [11] Hao Q, Chen G, Jeng MS. Frequency-dependent Monte Carlo simulations of phonon transport in two-dimensional porous silicon with aligned pores. J Appl Phys 2009;106:114321. [12] Tang GH, Bi C, Fu B. Thermal conduction in nano-porous silicon thin film. J Appl Phys 2013;114:184302. [13] Fu B, Tang GH, Bi C. Thermal conductivity in nanostructured materials and analysis of local angle between heat fluxes. J Appl Phys 2014;116:124310. [14] Chen G. Nanoscale Energy Transport and Conversion, A Parallel Treatment of Elections, Molecules, Phonons, and Photons. London: Oxford University Press; 2005. [15] He Y L, Xie T. Advances of thermal conductivity models of nanoscale silica aerogel insulation material. Appl Therm Eng, 2015;81:28-50. [16] Hsieh TY, Lin H, Hsieh TJ, et al. Thermal conductivity modeling of periodic porous silicon with aligned cylindrical pores. J Appl Phys, 2012;111:124329. [17] Pop E, Dutton RW, Goodson KE. Analytic band Monte Carlo model for electron transport in Si including acoustic and optical phonon dispersion. J Appl Phys 2004;96:4998-5005. [18] Asheghi M, Kurabayashi K, Kasnavi R, et al. Thermal conduction in doped single-crystal silicon films. J Appl Phys 2002;91:5079-5088. [19] Henry AS, Chen G. Spectral Phonon Transport Properties of Silicon Based on Molecular Dynamics Simulations and Lattice Dynamics. J Comput Theor Nanos 2008;5:141-152.

Biography Ya-Ling He is currently a professor of Xi’an Jiaotong University and director of the Key Laboratory of Thermo-Fluid Science & Engineering of MOE, China. She is an Associate Editor of Applied Thermal Engineering and Heat Transfer Research. She is also International Institute of Refrigeration-Vice President of Commission B1, and Members of the Scientific Council of ICHMT.