COMPUTATIONAL ASPECTS OF RADIATIVE TRANSFER EQUATION IN NON-ORTHOGONAL COORDINATES

Non-equilibrium energy transfer takes place for thin films when thermal disturbance is introduced. In this case, phonon transport inside the film governs the heat transport and temperature distribution in the film. In the present study an attempt is made to formulate and illustrate the phonon transfer in micro-scale silicon film of various shapes incorporating the non-orthogonal coordinate system. Successful application of the discrete-ordinates method to the solution of the equation for phonon radiative transport in non-orthogonal coordinates requires the application of various numerical techniques connected to the finite-difference method. The numerical solution of the equation for phonon transfer in non-orthogonal coordinate is introduced via adapting the discrete ordinate method. Phonon intensity distribution in the thin film is presented in terms of equivalent equilibrium temperature. It is found that film shape has significant effect on equivalent equilibrium temperature distribution inside the film. The validation study demonstrates that the code developed solving the equation for phonon transport is also applicable to the phonon transport in non-orthogonal coordinate system.


INTRODUCTION
Non-equilibrium energy transfer takes place when the size of the heat transferring becomes comparable to the mean free path of the solid material. In this case, the Fourier heat equation describing the heat transfer fails to predict correct temperature rise with in the small size solid [1]. The hyperbolic form of heat equation, such as Catteneo equation, does not fully describe the transport characteristics within the small size solid because of the wave nature of the energy transport. One of the models describing such energy transport incorporates the Boltzmann transport equation. After incorporating the phonon transport characteristics, particularly in small size dielectric films, the Boltzmann equation can be simplified and reduces to equation for phonon radiative transport (EPRT) [2]. Although the analytical solution of such equation is challenging, the attempt to have the analytical solution becomes viable for one-dimensional heating situation [3]. In the case of two-dimensional thin films, the analytical solution for the equation of phonon transport becomes extremely difficult because of the integrodifferential form of the governing equation. The numerical solution becomes viable incorporating the appropriate boundary conditions [4]. However, the general form of the equation can be applicable for all geometric configurations of the thin film, in general, the solutions are presented either rectangular or square films [5 -11]. This is because of the fact that the numerical solution for such geometry becomes easier to handle due to reflected phonons from the film edges. Consequently, further modifications are needed in the coordinate system to tackle such two-dimensional problems numerically.
The phenomenon of radiative transfer occurs in many physical situations such as the radiative heat transfer between surfaces through a participating media, conduction heat transfer through diffusive-ballistic phonons in thin films, neutron transport etc. In all these cases, the transport equations are very similar. In fact the left-hand-side of these equations is the same whereas the right-hand-side consists of terms specific to the phenomenon being described. So each equation is actually a radiative transfer equation (RTE). These equations can be solved analytically in closed form for simple situations. For problems involving complex geometries and other complications, the RTE may be solved numerically by various methods. One of the most common of such methods is the discrete-ordinates method (DOM). Until recently, the DOM method could only practically be applied to problems involving rectangular, cylindrical or spherical geometries. This was because the RTE was only available for such geometries. Vaillon et al. [12] derived the (RTE) for orthogonal coordinate systems but their derivation was incomplete. Building on the work of previous researchers, Freimanis [13] finally extended the radiative transfer equation to be applicable in any orthogonal or non-orthogonal coordinate system. Mansoor and Yilbas [14,15] employed this generalized RTE to the solution of phonon radiative transfer involving orthogonal coordinates.
In the author knowledge the RTE has not been solved by means of the DOM in non-orthogonal coordinates and the present study is an attempt to illustrate such an application by solving various problems related to the phonon transfer in micro-scale silicon film of various shapes. Successful application of the discrete-ordinates method to the solution of the equation for radiative transport in non-orthogonal coordinates, which requires the application of various numerical techniques connected to the finite-difference method.

MATHEMATICAL METHOD
The Equation of Phonon Radiative Transfer in non-orthogonal coordinates is written as [16], In Eq. (1), 1 2 3 ( , , ) x x x are the coordinate variables, t is the time variable, I  is the frequencydependent phonon intensity,  is the phonon frequency, v  is the phonon group velocity,   is the phonon mean-free-path,  is the polar angle and  is the azimuthal angle. The remaining variables 1 2 3 ,,    as well as the spatial derivatives of the polar and azimuthal angles are described below, In Equations (9) to (14), k ij  are the connection coefficients and are defined as [17], Finally, the derivatives of the polar and azimuthal angles are defined as [16],

APPLICATION TO STEADY-STATE, HEAT TRANSPORT IN MICRO-SCALE SILICON THIN FILMS
Consider a planar thin film of silicon of dimensions on the order of a micrometer or less. The flake may have an irregular, complex shape as shown in figure (1). For the solution of the Equation of Phonon Radiative Transfer (EPRT) this flake is now covered by a non-orthogonal grid. The grid can be generated numerically to obtain the following relations, The physical domain is described by the xz − coordinate system and the computational domain is described by the  − coordinate system as shown in figure (1).

Figure 1.
Coordinate system and grid used in orthogonal and non-orthogonal arrangements.
The two-dimensional, steady-state, frequency-independent EPRT for the above described situation is obtained from Eq. (1) by making the following assignments, 13 , xx  →→ , With the following expressions for the coefficients, Eq. (21) may now be solved numerically by approximating its various derivatives by finite difference formulas as well as by approximating the integral on the right-hand-side. Much of the details of the numerical method can be obtained from [14] and [15]. Here, we discuss some details that are specific to the solution of the EPRT in a non-orthogonal coordinate system. The first to note is that in general the sign of the coefficients of all the derivatives is a priori unknown. Hence, to obtain a stable discretization, streamwise diffusion has to be introduced, i.e. an upwind differencing scheme has to be used, in general, for all derivatives. This means that for any derivative in Eq. (20), we use a backward difference formula at a point where its coefficient is positive and a forward difference formula where its coefficient is negative. The second point to note is that the upwind discretization works well for the derivatives I   Figure (2) shows the two-dimensional silicon films of four different shapes. The non-orthogonal grids incorporated in the numerical simulations covering are also shown for each film. In order to quantify the phonon intensity distribution inside the film, equivalent equilibrium temperature is introduced. Equivalent equilibrium temperature is associated with the average energy of all phonons around a local point when the phonon energies are redistributed adiabatically to an equilibrium state. In addition, in the simulation, the steady-state, frequencyindependent equation for phonon transport is solved when thermal disturbance is introduced from the edges of 167 the film. Hence, the thermal disturbance is introduced through the constant temperature boundary conditions, i.e. the bottom boundary is maintained at a temperature of 301 K whereas all the remaining boundaries are maintained at a temperature of 300 K.

Figure 2.
The silicon film with different shapes are shown together with non-orthogonal grids used in the simulations Figure (3) shows the contours of the equivalent equilibrium temperature for the four different silicon film geometries. High temperature edge (301 K) of the film remains at the film bottom and other edges are set to be 300 K. Temperature distribution varies significantly in each film. This indicates that the shape factor for small scale films is extremely important in terms of heat transfer; in which case, temperature extends in the film bottom particularly for "M" shape film. This is because of the fact that the phonons emanating from the high temperature edge undergoes scattering prior to reaching at the low temperature edge. The phonons, which are reflected from the edges, also undergo scattering and contributing to the equivalent equilibrium temperature rise in the film. However, low temperature edges act like a phonon sink while reducing the phonon intensity in the near region of the low temperature edges. Therefore, equivalent equilibrium temperature remains low in the region of the cold edges. The low temperature extends considerably in the film in the case of trapezium and "M" shape films. In order to validate and compare the non-orthogonal grid formulation, two film shapes are considered, namely square and quadrilateral, as shown in figure (4). The boundary conditions are set to be same as those shown in figure (3). However, orthogonal grids are used in the solution of the equation for phonon transport. Hence, a numerical code is developed for this purpose. Figure (4) shows the equivalent equilibrium temperature contours for the square shape film as well as that of the quadrilateral shape film. Again, for the square shaped film, the equation for phonon transport is solved numerically via a code developed for rectangular geometries only. However, for the quadrilateral flake, the equation for phonon transport is solved numerically via a code developed for irregular shapes but covered with the orthogonal grid arrangements. It is evident that the equivalent equilibrium temperature contours are same as that obtained from the equation for phonon transport code for the non-orthogonal grids as shown in figure (3). This similarity validates not only the code developed to solve the equation for phonon transport, but it validates the equation of phonon transport that is also applicable to non-orthogonal coordinate system.

CONCLUDING REMARKS
Since the equation for phonon transport in non-orthogonal coordinates has not been solved numerically by means of the discrete ordinate method. Consequently, in the present study, the formulation of equation of phonon transport in non-orthogonal coordinate system is introduced. An attempt is made to illustrate such application of equation for phonon transport via solving various problems related to the phonon transfer in the micro-scale two-dimensional silicon film of various shapes. It is demonstrated the successful application of the discrete-ordinates method to the solution of the equation for radiative phonon transport in non-orthogonal coordinates. This requires the application of various numerical techniques connected to the finite-difference method. Moreover, film shapes have significant effect on the equilibrium equivalent temperature distribution inside the film. The phonons emitted from high temperature edge undergo scattering in the film while causing equivalent equilibrium temperature rise in the film. The low temperature edges act like phonon sinks while lowering equivalent equilibrium temperature in the film. Therefore, low temperature region extends in the region close to these edges, which is more pronounced for "M" shape. The validation study demonstrates that the code developed solving the equation for phonon transport is also applicable to the phonon transport in non-orthogonal coordinate system.

ACKNOWLEDGMENTS
The author acknowledges the support of King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia, for this work.