Hybrid discrete ordinates—spherical harmonics solution to the Boltzmann Transport Equation for phonons for non-equilibrium heat conduction

https://doi.org/10.1016/j.jcp.2011.05.024Get rights and content

Abstract

The Boltzmann Transport Equation (BTE) for phonons has found prolific use for the prediction of non-equilibrium heat conduction phenomena in semiconductor materials. This article presents a new hybrid formulation and associated numerical procedures for solution of the BTE for phonons. In this formulation, the phonon intensity is first split into two components: ballistic and diffusive. The governing equation for the ballistic component is solved using two different established methods that are appropriate for use in complex geometries, namely the discrete ordinates method (DOM), and the control angle discrete ordinates method (CADOM). The diffusive component, on the other hand, is determined by invoking the first-order spherical harmonics (or P1) approximation, which results in a Helmholtz equation with Robin boundary conditions. Both governing equations, referred to commonly as the ballistic-diffusive equations (BDE), are solved using the unstructured finite-volume procedure. Results of the hybrid method are compared against benchmark Monte Carlo results, as well as solutions of the BTE using standalone DOM and CADOM for two two-dimensional transient heat conduction problems at various Knudsen numbers. Subsequently, the method is explored for a large-scale three-dimensional geometry in order to assess convergence and computational cost. It is found that the proposed hybrid method is accurate at all Knudsen numbers. From an efficiency standpoint, the hybrid method is found to be superior to direct solution of the BTE both for steady state as well as for unsteady non-equilibrium heat conduction calculations with the computational gains increasing with increase in problem size.

Introduction

The efficient removal of heat from modern-day solid-state electronic and optoelectronic devices is a daunting task, and overheating is one of the most common causes of device failure. The mean free path of the energy-carrying acoustic wave packets (or phonons) in silicon at room temperature is approximately 300 nm [1]. On the other hand, characteristic dimensions of modern semiconductor devices range from a few tens of nanometers to a few hundreds of nanometers [2]. Consequently, heat conduction in such devices cannot be described adequately using continuum equations, namely the Fourier law of heat conduction. The Boltzmann Transport Equation (BTE) for phonons has been used successfully in the past for the prediction of non-equilibrium heat conduction phenomena, and continues to be popular because of its validity over a large range of length scales [2]. Simulation of a practical electronic device requires coupling of charge carrier transport with phonon transport in order to address electron-phonon coupling and related phenomena. Although significant progress has been made in the simulation of charge carrier transport in semiconducting materials, comparable progress has not been made on the thermal side, i.e., in modeling phonon transport. The present article contributes to this longstanding need.

In general, the BTE is a seven-dimensional nonlinear integro-differential equation: 3 spatial coordinates, 3 wave vector coordinates, and time. Under the so-called single time relaxation approximation, it reduces to a seven-dimensional linear partial differential equation. Even if the physics of phonons were to be neglected, i.e., no dispersion, no polarization, frequency-independent (or gray), the BTE is a six-dimensional partial differential equation: 3 spatial coordinates, 2 angular coordinates, and time. This high dimensionality of the BTE makes it very challenging to solve. Traditionally, stochastic methods, such as the Monte Carlo method, have been effectively employed for the solution of such high-dimensional partial differential equations. For example, the Monte Carlo method has been successfully used to solve the BTE or its variants for rarefied gas (also known as direct simulation Monte Carlo or DSMC) [3], photons (or radiative transport) [4], [5], [6], electrons and other charge carriers [7], [8], [9], and the stochastic differential equations resulting from modeling turbulent flow and related phenomena [10], [11]. The advantage of the Monte Carlo method is that it is tractable and almost linearly scalable for large number of dimensions, and is amenable to addressing complex physics via interactions between the stochastic samples. The shortcoming of this method is that it is prohibitively expensive for practical engineering applications. In the case of phonon transport, for example, studies [12], [13], [14] have shown that even simulation of heat conduction in a one-dimensional thin film can require several tens of hours of CPU time. In light of this shortcoming, it is fair to contend that the Monte Carlo method is useful for generating benchmark solutions for simple problems but not practical for simulating non-equilibrium heat conduction in large-scale three-dimensional structures, and deterministic methods for solving the BTE are desired.

Deterministic solution of the BTE has its roots in the neutron and photon (radiation) transport area. The Discrete Ordinates Method (DOM) [15] based on the SN approximation has been used extensively for the computation of photon (radiation) transport. It has been adopted for phonon transport and brought to the limelight primarily by Murthy and co-workers [16], [17], [18], [19]. Over the years, the core DOM has been refined to improve accuracy and efficiency. For example, methods have been developed to enable specular and partially specular reflection at boundaries [16], [20]. Also, in an effort to reduce so-called “ray effects” and “false scattering”—two well-known shortcomings of the discrete ordinates method—a switch has been made from the standard DOM to the so-called Control Angle Discrete Ordinates Method (CADOM) [21], [22], [23]. Higher-order and more robust discretization schemes in space, such as the SMART scheme [24], has been adopted from the computational fluid dynamics literature to improve accuracy and convergence [17], [18] of the phonon BTE. While the discrete ordinates method and its variants have shown tremendous promise for solution of the BTE for phonons, it is computationally quite expensive. In particular, in the ballistic regime (high Knudsen number), a large number of directions (or control angles) have to be used to attain acceptable accuracy, as has been shown in previously published results [17], [18], [19], and will also be shown in this article. This implies that solution of a large number of partial differential equations is necessary. On the other hand, in the diffusive regime (low Knudsen number), since phonon transport is diffusive, use of DOM with high angular resolution is wasteful and use of a simpler diffusive approximation is warranted for improved efficiency.

One method which has often been used to solve the BTE is the method of spherical harmonics (or PN approximation). In this method, the angular directions are not discretized. Instead, spherical harmonic basis functions, namely Legendre polynomials, are used to capture angular variations in the intensity through analytical integration. The lowest order spherical harmonics approximation, namely the P1 approximation, reduces the BTE to a single Helmholtz equation with Robin boundary conditions, making it an attractive choice for solution of the BTE since only one partial differential equation has to be solved, as opposed to several tens in the discrete ordinates method. Unfortunately, the P1 approximation is accurate only in the diffusive (low Knudsen number) regime, and its accuracy at intermediate or high Knudsen numbers is unacceptable, as conclusively demonstrated for photon (radiation) transport [25], [26]. Therefore, the P1 approximation is inappropriate for solution of the BTE in truly non-equilibrium scenarios (intermediate to high Knudsen numbers) where the angular variation of the intensity is significant. Higher order PN approximations, as recently formulated for radiation transport [26], [27], [28], may be used to alleviate this problem. However, use of higher order PN approximations is significantly more cumbersome than using the P1 approximation [26], and the benefits in terms of improved accuracy is often marginal [28].

The so-called Modified Differential Approximation (MDA) was proposed to remove the shortcomings of the P1 approximation for intermediate and high Knudsen numbers. In this method, first proposed by Olfe [29], and later generalized for radiation transport by Modest [30], the intensity of the energy carrier is split into two components: ballistic and diffusive. The ballistic component is determined using a surface-to-surface exchange formulation that employs geometric viewfactors, while the spherical harmonics approximation is invoked for the diffusive component, for which it is justifiably applicable. The result is a hybrid approach that is expected to be efficient and accurate at all Knudsen numbers. In the past decade, the MDA approach has been adopted for solution of the BTE for phonons by Chen and co-workers [31], [32], [33], resulting in the so-called ballistic-diffusive equations (BDE) of phonon transport. In the BDE formulation proposed by Chen and co-workers [31], [32], [33], the specific heat capacity of the material appears as an input. This makes direct comparison with results of the BTE difficult, since the specific heat capacity of the material is not an input in the BTE. Secondly, Chen and co-workers introduce artificial temperatures, namely “ballistic” and “media” temperatures, in their formulation. These temperatures do not have physical meaning and are introduced as mathematical artifacts. As a result, they make the formulation—in particular, the boundary conditions—difficult to understand and interpret. A final point to note is that the surface-to-surface exchange formulation used by Chen and co-workers [31], [32], [33], which employs geometric viewfactors, is prohibitive on two counts. First, it is very expensive and tedious for complex multi-dimensional geometries in which case determination of the viewfactors itself is a monumental task [34]. Secondly, all surfaces in a heat conduction simulation are not necessarily diffuse, and therefore, the use of diffuse geometric viewfactors is limiting in its scope.

In this paper, we employ a recently developed [35], [36] new hybrid formulation (or BDE formulation) for approximate solution of the BTE. The new formulation does not require a priori knowledge of the specific heat capacity of the material, and requires the exact same inputs as the original BTE, making direct comparison with the BTE possible. While the new formulation is based upon the same intensity splitting philosophy originally proposed by Olfe [29] for development of the MDA formulation, the resulting BDE is solved using procedures that are general enough to be applicable to any arbitrary geometry with arbitrary thermal boundary conditions. The ballistic component is determined, in this case, using the Discrete Ordinates Method (DOM) and its variants such as the Control Angle Discrete Ordinates Method (CADOM), while the diffusive component is determined by invoking the P1 approximation. The resulting governing equations and boundary conditions are discretized on an unstructured mesh using the finite-volume procedure. The final set of algebraic equations, both for the ballistic and diffusive components, is solved using an incomplete LU (ILU) pre-conditioned GMRES solver [37]. As part of this study, direct comparisons are made with benchmark results obtained using the Monte Carlo method, as well as standalone DOM or CADOM for a benchmark transient two-dimensional heat conduction problem at various Knudsen numbers. Subsequently, the code is explored for complex two-dimensional (2D) and three-dimensional (3D) geometries. In addition to investigation of accuracy, convergence and timing studies are conducted for various mesh sizes and angular resolutions. While the results shown here are for the gray BTE, the formulation is applicable to the non-gray (frequency dependent) BTE with minor modifications. While Ref. [36] presents the mathematical formulation of the proposed hybrid approach in great detail and presents only simple 2D test cases as preliminary proof of concept, the present paper is aimed toward development of the numerical methods required to solve the resulting governing equations in complex 3D geometries, and elucidating the numerical issues, namely convergence, efficiency, and memory requirements associated with solution of the governing equations.

Section snippets

Theory and mathematical formulation

Quantized lattice vibrations or phonons are the predominant carriers of thermal energy in semiconductor materials [38]. If the mean free path of the traveling phonons is larger than the characteristic dimension of the device being modeled, thermodynamic equilibrium ceases to exist, and thus, the Fourier law of heat conduction is invalid. Heat conduction, in such a scenario, is referred to as non-equilibrium heat conduction.

The Boltzmann Transport Equation (BTE) is a semi-classical equation, and

Numerical procedure

The preceding sections outline the governing equations that need to be solved and the associated boundary and initial conditions for the BDE. In order to attain the solution to a non-equilibrium heat conduction problem using this formulation, the equations need to be discretized and solved. The procedures adopted to do so are described below.

Results and discussion

In order to test the new hybrid formulation and algorithms described in the preceding sections, three different unsteady heat conduction problems were considered. These test cases were chosen with the objective of assessing both accuracy and efficiency of the method. The geometry of the first test case was chosen to be a two-dimensional square cavity so that benchmark Monte Carlo results, which are extremely time-consuming to compute, could be generated with relative ease. On the other hand,

Summary and conclusions

A new generalized form of the Ballistic-Diffusive Equations (BDE) for approximate solution of the Boltzmann Transport Equation (BTE) for phonons is demonstrated. This new formulation, unlike previously published formulations of the BDE [31], [32], does not require a priori knowledge of the specific heat capacity of the material. The only input required in this formulation is the phonon scattering time scale (or Knudsen number in non-dimensional form), which is the same input required for

Acknowledgments

Financial support for this work was provided in part by the Department of Energy’s Basic Energy Science Program through Grant Number DE-FG02-06ER46330. ESI Group is acknowledged for providing free licenses to their mesh generation package, CFD-GEOM™.

References (43)

  • S. Mazumder et al.

    A fast Monte-Carlo scheme for thermal radiation in semiconductor processing applications

    Numer. Heat Transfer, Part B: Fundam.

    (2000)
  • J.D. Maltby et al.

    Performance, accuracy, and convergence in a three-dimensional Monte Carlo radiative heat transfer simulation

    Numer. Heat Transfer, Part B: Fundam.

    (1991)
  • M.V. Fischetti et al.

    Monte Carlo study of electron transport in silicon inversion layers

    Phys. Rev. B

    (1993)
  • C. Jacoboni et al.

    The Monte Carlo method for the solution of charge transport in semiconductors with applications to covalent materials

    Rev. Mod. Phys.

    (1983)
  • S. Mazumder et al.

    Monte Carlo study of phonon transport in solid thin films including dispersion and polarization

    J. Heat Transfer

    (2001)
  • D. Lacroix et al.

    Monte Carlo transient phonon transport in silicon and germanium at nanoscale

    Phys. Rev. B

    (2005)
  • A. Mittal et al.

    Monte Carlo study of phonon heat conduction in silicon thin films including contributions of optical phonons

    J. Heat Transfer

    (2010)
  • W. Fiveland

    Three-dimensional radiative heat transfer solutions by the discrete ordinates method

    J. Thermophys. Heat Transfer

    (1988)
  • J.Y. Murthy et al.

    Computation of sub-micron thermal transport using an unstructured finite-volume method

    J. Heat Transfer

    (2002)
  • S.V.J. Narumanchi, Simulation of Heat Transport in Sub-Micron Conduction, PhD. Dissertation, Department of Mechanical...
  • S.V.J. Narumanchi et al.

    Sub-micron heat transport model in silicon accounting for phonon dispersion and polarization

    J. Heat Transfer

    (2004)
  • Cited by (36)

    • Improvements to a class of hybrid methods for radiation transport: Nyström reconstruction and defect correction methods

      2020, Journal of Computational Physics
      Citation Excerpt :

      One common decomposition approach involves splitting the radiation distribution function into a sum of two components – one containing non-collisional particles and one containing collisional particles – both of which exist over the entire problem domain. Methods based on this approach have been used in a variety of contexts under different names, including first-collision source methods [6,7], ballistic-diffusive models [8–15], and the isotropic diffusion source approximation (IDSA) [16–20]. While the contexts under which these methods are applied may differ, the same approach is used in each case: a method tailored toward the free-streaming limit, such as a ray-tracing approximation, is applied to the non-collisional component, and a method tailored toward the strongly-collisional limit, such as a diffusion approximation, is applied to the collisional component.

    • Fractional parabolic two-step model and its accurate numerical scheme for nanoscale heat conduction

      2020, Journal of Computational and Applied Mathematics
      Citation Excerpt :

      But it is not straightforward to solve the BTE in general because of its high dimensionality [22] . Although many methods including the Monte Carlo method [27–31], the lattice Boltzmann method [32,33], and various deterministic discretization-based methods [34–42] have been developed to numerically solve the BTE, it is not straightforward to include scalar laser field in MD or phonon BTE. Therefore, it is of great interest to develop innovative and accurate models with less computational cost and more simplicity for analyzing nanoscale heat transfer induced by ultrafast pulsed laser heating.

    View all citing articles on Scopus
    View full text