Flexible and Accurate Simulation of Radiation Cooling with FETD Method

Thermal management and simulation are becoming increasingly important in many areas of engineering applications. There are three cooling routes for thermal management, namely thermal conduction, thermal convection and thermal radiation, among which the first two approaches have been widely studied and applied, while the radiation cooling has not yet attracted much attention in terrestrial environment because it usually contributes less to the total amount of thermal dissipation. Thus the simulation method for radiation cooling was also seldom noticed. The traditional way to simulate the radiation cooling is to solve the thermal conduction equation with an approximate radiation boundary condition, which neglects the wavelength and angular dependence of the emissivity of the object surface. In this paper, we combine the heat conduction equation with a rigorous radiation boundary condition discretized by the finite-element time-domain method to simulate the radiation cooling accurately and flexibly. Numerical results are given to demonstrate the accuracy, flexibilities and potential applications of the proposed method. The proposed numerical model can provide a powerful tool to gain deep physical insight and optimize the physical design of radiation cooling.

adopted numerical methods for the transient thermal analysis include integral equation-based method, finite-difference time-domain (FDTD) method and finite-element time-domain (FETD) method [10][11][12][13][14][15][16][17][18][19][20][21][22] , etc. The integral equation-based method is very suitable for open problems due to the involvement of Green's function. But its level of complexity will increase enormously when dealing with inhomogeneous material. The differential equation-based methods, such as FDTD and FETD, are preferred choice for analyzing complex inhomogeneous material in bounded geometries. The FDTD method is widely applied due to its simplicity and computational efficiency. Since orthogonal grids are employed in the conventional FDTD method, leading to the staircase error. Consequently, the mesh density must be increased dramatically to model objects with curved surfaces. By contrast, the FETD method with tetrahedral elements is more flexible to model arbitrarily-shaped geometries. Most of the aforementioned studies on numerical methods focus on the simulation of conduction cooling and convection cooling. Very few of them involve the simulation of radiation cooling. In particular, the radiation boundary condition used in these studies is an approximate one which neglects the wavelength and angular dependence of the emissivity of the object surface. For the real-world applications, this emissivity is always a function of the wavelength and direction angle. In this paper, the heat conduction equation with a rigorous radiation boundary condition is discretized with the FETD method to simulate the radiation cooling phenomena accurately and flexibly. Fig. 1, consider an arbitrarily shaped object Ω within which there is no bulk motion, the governing equation for its transient thermal analysis can be expressed as

Governing Equation and Its Definite Conditions. As shown in
where ρ is the density of the materials, c ρ is the specific heat capacity, T denotes the temperature distribution as a function of time and space, κ is the thermal conductivity, and Q represents the heat source, corresponding to the energy generated per unit volume in the medium. The differential Eq. (1) with predefined boundary condition and initial condition can be solved to obtain the spatial temperature distribution in the medium. There are four categories of commonly-used boundary conditions in thermal analyses. For the convenience of description, the surface of Ω is denoted by S = S 0 ∪ S 1 ∪ S 2 ∪ S 3 . If S 0 is maintained at a fixed temperature T s , the Dirichlet boundary condition will be applied: s 0

= ∈
The second category of boundary condition is termed as the Neumann condition, corresponding to a fixed heat flux at the surface S 1 : where n denotes the unit outward normal vector of the boundary surface. The surface S 1 becomes a perfectly insulated or adiabatic surface when q s = 0. The third category of boundary condition is the convection boundary condition obtained from the Newton's law of cooling: where h is the convective heat transfer coefficient. T sur refers to the surrounding temperature. This boundary condition is utilized to describe the heat transfer between the surface S 2 and the fluid moving over the surface. If the temperature of the medium is above absolute zero, it will radiatively emit heat to outer space, which is called thermal radiation. The fourth category of boundary condition, namely radiation boundary condition, will be applied if thermal radiation occurs between the surface S 3 of the medium and its surrounding environment: where P rad (T) denotes the power radiated by the structure per unit area. P atm (T sur ) represents the power absorbed from the surrounding atmosphere per unit area. ε(λ, θ, ϕ) is the spectral and directional emissivity of S 3 . As shown in Fig. 2, θ is the zenith angle, while ϕ is the azimuthal angle. λ denotes the wavelength. I λ,B (λ, T) refers to the spectral radiance of a blackbody at temperature T.
( ) where c 0 = 2.998 × 10 8 m/s is the speed of light in free space. k B = 1.381 × 10 −23 J/K is the Boltzmann constant. h = 6.626 × 10 −34 J · K is the Planck constant. Equations (6) and (7) are the rigorous expression for P rad (T) and P atm (T sur ). In the traditional radiation boundary condition, the emissivity of the medium surface is assumed to be independent of the wavelength and direction, namely ε(λ, θ, ϕ) = ε 0 . Since the total emissive power of a blackbody can be calculated as where σ = 5.68 × 10 −8 W/m 2 K 4 is the Stefan-Boltzmann constant. We can obtain the traditional radiation boundary condition as When adopting the rigorous expression of radiation boundary condition in (5), the trapezoidal rule of numerical integration is utilized to calculate (6) an (7) 23 . It is worth mentioning that the infinite integral over the wavelength should be truncated into a finite one. Equation (9) can be used to test the accuracy of different truncation boundary and resolution. Numerically, we find that 0.4% relative accuracy can be achieved with a spectral resolution of 2 nm from 1 μm to 1000 μm and with 3 degree angular resolution.
In Cartesian coordinates, the initial condition can be specified as Finite-Element Time-Domain Method. The FETD method is employed to solve Eq. (1) with the above mentioned boundary conditions and the initial condition for obtaining the time-varying temperature distribution of the medium. Firstly, the medium Ω is discretized into a set of tetrahedral elements Ω i . Nodal basis functions are employed to expand T in each element Ω i

24
. Then considering the aforementioned Neumann, convection and radiation boundary conditions simultaneously, we can obtain the spatially discrete form of equation (1) by using the Galerkin's scheme: N i and N j refer to the i-th and j-th nodal basis functions, respectively. {T} refers to the time-dependent expansion coefficients to be solved, which is also the vector constituted by the temperature at the nodes of the elements. For the nodes residing on S 0 , we need to impose the Dirichlet boundary condition. The corresponding unknown expansion coefficients on S 0 are explicitly given by (2), only the unknown coefficients not over S 0 are to be solved by (12). Suppose that the temperature of the nodes at time t i = iΔt is {T i }. Δt is the time step size. The Crank-Nicolson scheme is adopted for the temporal discretization of (12) to arrive at an unconditionally stable system 25 : Note that (16) is a nonlinear system of equations. Particularly, the last term arising from the radiation boundary condition (5) requires to compute the radiation power P rad (T) at the time t = t i . In order to avoid solving nonlinear equation, Eq. (16) at the time t = t i is solved in an iterative manner: (17), which provides the values for the right hand side of (17), leaving the temperature vector ( 1) as unknowns. Then ( 1) can be solved by using the multifrontal method for sparse linear equations 26 . This process is repeated until the difference between ( ) is less than an acceptable tolerance. It is set to be 10 −5 in this paper. Finally, we can obtain = 1) . Figure 3 shows the schematic picture for the numerical solution of the temperature distribution at each time step.

Results
Silicon Cube. As shown in Fig. 4(a), a silicon cube with the side length of 0.5 m is studied to validate the accu- . The ambient temperature is set to be 300 K. The initial temperature of the cube is 800 K. The cube is meshed into 3736 tetrahedral elements as shown in Fig. 4(b). The temperature at a randomly chosen observation point r = (0.185 m, 0.18 m, 0.256 m) is recorded in the following simulation.
Firstly, the Dirichlet, Neumann, convection and traditional radiation boundary conditions are imposed to the whole boundary surface, respectively. The temperatures on the boundary surface are fixed to be 300 K in the Dirichlet boundary condition. The fixed heat flux in Neumann boundary condition equals to 800 K/m 2 . The convective heat transfer coefficient is set to be h = 15 W/[m 2 · K] in the convection boundary condition. The emissivity of the boundary surface is set to be ε 0 = 0.9 in the traditional radiation boundary condition. Both the proposed method and the COMSOL software are employed to analyze this problem 27 . The temporal temperature at the observation point is shown in Fig. 5(a). Very good agreement between the proposed method and the COMSOL software is achieved. We can also find that the temperature at the observation point before 14000 s decreases faster Scientific REpoRTS | (2018) 8:2652 | DOI:10.1038/s41598-018-21020-w when applying the radiation boundary conditions than using the convection boundary condition, which implies that the radiation cooling can be very effective especially when the temperature of object surface is high since it is proportional to the fourth power of temperature.
Secondly, different boundary conditions are imposed to different surfaces of the cube to further check the code. The fixed outward heat flux 2000 W/[m 2 ] is applied on the surfaces of x = 0 and x = 0.5. The convection boundary condition is imposed on the surfaces of y = 0 and y = 0.5. The radiation boundary condition is applied on the surfaces of z = 0 and z = 0.5. As shown in Fig. 5(b), four groups of the convective heat transfer coefficients and emissivities are tested, again good agreement between COMSOL and the proposed method can be observed. Moreover, it is easy to find that a larger convective heat transfer coefficient or emissivity can lead to faster cooling.
Thus far, we have validated the implementation of the Dirichlet, Neumann, convection and traditional radiation boundary conditions. However, the treatment of the rigorous radiation boundary condition in (5) has not been proved yet. According to (9), if we set the emissivity to be a constant ε 0 independent of the wavelength and direction, the results obtained through utilizing (5) and (10) should be the same. This can be used to check the correctness of treatment of the rigorous radiation boundary condition. Figure 5(c) shows the results  corresponding to ε 0 = 0.2, ε 0 = 0.5 and ε 0 = 0.9 when applying the radiation boundary condition to the six surfaces of the cube. Obviously, the results of the rigorous radiation boundary condition agree well with those of the traditional radiation boundary condition, which demonstrates the correctness of our proposed method.
Actually, the emissivity of an object surface is seldom uniform, and it is usually dependent on the wavelength and emission angle. In order to exhibit such a feature, we consider three different kinds of emissivities: (1) The emissivity only depends on the wavelength. It equals to 0.9 in the wavelength range from 8 μm to 11 μm and is zero elsewhere. (2) The emissivity only depends on the emission angle. It equals to 0.9 in the angle range from 0 to π/3 and is zero in the angle range from π/3 to π/2. (3) The emissivity depends both on the wavelength and emission angle. It equals to 0.9 in the wavelength range from 8 μm to 11 μm and in the angle range from 0 to π/3. The emissivity is zero elsewhere. The rigorous radiation boundary condition is applied with the above three different kinds of emissivities. The temperature of the observation point can be obtained by the proposed method as shown in Fig. 5(d). Note that the case of ε 0 = 0.9 in Fig. 5(c) is also replotted in Fig. 5(d) for comparison. Obviously, the cooling effect will be overestimated if adopting the traditional radiation boundary condition. It is worth mentioning that the number of iterations for the solution of (17) at each time step in these three cases are less than 5.
3D IC Package with a Heat Sink. Consider a 3D IC package model with a heat sink, the cross-section view of which is shown in Fig. 6(a). The model is composed of a heat sink, thermal interface materials (TIM), a chip, 6 × 6 through-silicon-vias (TSVs) and an underfill layer. Its geometry and material parameters are summarized in Table 1. A uniform volume power of 1 W is assumed on the chip layer as the heat source. Initially, the system temperature is supposed to be room temperature of 300 K. To perform the numerical simulation, this model is discretized into 269,446 tetrahedral elements. The temperature of a randomly chosen observation point r = (7.25 mm, 0.7 mm, 2.875 mm) at the top surface of the chip is recorded.
Case 1: Firstly, the convection boundary condition is imposed on the surface of the heat sink with h = 15W/ [m 2 · K] to model the natural convection. Other surfaces are assumed to be adiabatic. Both the proposed method Case 2: Actually, anodic oxidation treatment is usually applied to the surface of the heat sink to improve its emissivity. Thus the effect of radiative cooling can be improved. Suppose that the surface emissivity of the heat sink equals to ε 0 = 0.93 after the anodic oxidation treatment. Natural convection with h = 15 W/[m 2 · K] still exists. Then the heat conduction equation combined with the radiation boundary condition and convection boundary condition can be utilized to analyze this problem. Figures 6(b) and 7(c,d) show that the results obtained by the proposed method agree well with those by COMSOL software. The highest temperature of the observation point decreases nearly 20 °C with the radiative cooling.
Case 3: In real-world applications, the assumption of a wavelength and direction independent emissivity is often invalid. The surface emissivity usually strongly depends on the wavelength and direction, which can be measured and then used for accurate simulation of radiation cooling. Suppose that the surface of the heat sink is treated to have a measured emissivity as shown in Fig. 6(c), corresponding to an average emissivity of 0.38. The cases of supposed emissivity and average emissivity are simulated based on the rigorous radiation boundary condition and the traditional radiation boundary condition, respectively. As shown in Figs 6(b) and 7(e,f), the result of the traditional wavelength and direction independent boundary condition apparently deviates from the   28 . The metamaterial is a 50 μm thick film containing 6% of microspheres by volume which has an averaged infrared emissivity 9.3. This material has a promising prospect for passive radiative cooling. Suppose the heat sink in this  example is wrapped with this film and thus its surface emissivity will be improved. The measured emissivity of the film in reference 28 is replotted in Fig. 6(c). Then the temperature distribution of the 3D IC package can be analyzed again by our proposed method. The temporal temperature of the observation point is given in Fig. 6(b). Since the measured emissivity of the film has a weak dependence on the wavelength and is very closed to the ideal wavelength and direction independent emissivity, the cooling effect with the practical emissivity is also very close to that with the ideal wavelength and direction independent emissivity. For the aforementioned four cases, the relative errors of the temperature at the observation point between the proposed method and COMSOL software are shown in Fig. 6(d). The CPU time for the whole simulation, peak memory and average number of iterations per time step are given in Table 2. All of the numerical experiments are performed on 2.4 GHz CPU and 8 GB RAM.