# Modelling and Numerical Simulations of Heat Distribution for LED Heat Sink

Authors: Jianxin Zhu; Lixia Sun
Journal: Discrete Dynamics in Nature and Society
 (2016)

---

## Abstract

Light-emitting diode (LED) has higher efficiency and longer lifetime when compared with the conventional lighting. However, the efficiency and lifetime will be degraded greatly when it is operated at a high temperature. Now, both previous simulation and experimental results have already indicated that the heat transfer in vertical direction of the LED lamp by conduction is the most critical component. In this paper, a simplified numerical simulation model is built to estimate the heat distribution of the LED heat sink in the spherical coordinate system, which would be useful for its shape optimization design. With this model, some mathematical treatments are provided to a heat conduction equation, in order to rapidly compute the static heat distribution and the temperature of different designs of LED heat sinks. The built rapid heat sink evaluation method, implicit finite difference method (IFDM), is unconditionally stable. Several heat distribution simulations could demonstrate that our built mathematical model conforms well to the reality and our method is full of feasibility and effectiveness.


---

## Body

## 1. Introduction

Recently, light-emitting diode (LED) lamp has played an important role in the illumination market, mainly due to the advantages of more compact package in size, higher efficiency, and longer lifetime than the conventional lighting types. It has become the fourth-generation light source in the auto industry [1, 2]. Moreover, it has been considered as the ultimate light source. However, when we want to produce more light output, the higher power needs to be applied to LED lamp. Subsequently, the amount of heat generated from the LED lamp will be greatly increased at the same time, which could reduce both the lifespan and the luminous efficiency. Researchers have found that both the efficiency and lifetime will be degraded pronouncedly when the LED is operated at a temperature higher than 70∘C. Thus, the design of heat sink used to cool LED lamp is full of importance to ensure their long lifetime and high efficiency.Numerous heat sinks in different shapes, sizes, and structures are currently used in the commercial LED light bulbs. However, most of these heat sinks were not optimized in cooling efficiency, dimension, and material cost. Lots of researchers have already studied the cooling of LED lamps [3–20], which are mainly operated under natural or forced convection. More specifically, a numerical and experimental research considering the whole set constituted by the vapour chamber and a finned heat sink was conducted in [3]. A 3D finite element simulation was presented for an array of high power LED lamp with a heat sink in [4], where the thermal resistance network was performed to estimate different contributions for the heat management. In [5], a radial heat sink with pin fins was optimized while keeping a similar cooling performance as in the previous study. In the related further work [17], a three-dimensional fin-height profile heat sink was optimized. In [11], the enhanced model with an aluminum nitride (AIN) insulation plate was more effectively conducted than the conventional chip on board (COB) model. In [12], LED effectively produced a thermal effect by the changes of the fin shape and intervals of various heat sink designs. In [13], the thermal performance of LED illumination systems could be improved by a detailed CFD analysis. In [14], the numerical simulation and optimization of a radial heat sink were performed, and the impossibility of optimizing both thermal performance and heat sink mass at the same time was concluded. In [15], some experimental results were shown for a similar heat sink configuration. In [16], the effect of both radiation and natural convection on the thermal performance of a radial heat sink was studied similar to [14, 15]. In [20], numerical simulation was conducted to investigate the impingement and film composite cooling on blade leading region. The luminance of LED lamp is greatly related to junction temperature, so evaluating the junction temperature is one of the methods to determine its heat dissipation [6–10]. In [6], the junction temperature was assessed by the heat dissipation simulation of a single chip high power LED package, and an electric-heat-optical system dynamics model for LED luminance control was provided in [7]. A 3 W high power LED array system with an in-line pin fin heat sink was presented in [8]. In [9], a tandem 12-chip module with three types of heat sinks was discussed, while the impact factors of the junction temperature were studied in [10] through the thermal performance of conventional plate-fin heat sinks and novel cooling device of heat conductive plates. Recently, the related research can also be found in [18, 19].From the above research, we could find that it is major importance to get the heat distribution when we want to optimize the designs of LED heat sinks. To the LED heat sink, recent simulation and experimental results have already proved that the heat transfer in vertical direction of the LED lamp by conduction is the most critical component. Then we mainly consider the heat conduction equation concerning the LED heat sink in this paper. Many methods have been developed to the heat transfer problems in the Cartesian coordinate system, such as the technique of alternating directional implicit method (ADIM) [21], the finite volume method (FVM) [22], and the finite element method [23]. However, for the LED heat sink, its model of heat distribution can be much simpler in the spherical coordinate system than in the Cartesian coordinate system, though there is a little research about 3D heat transfer problems in the spherical coordinate system. The objective of our present work is building a simplified mathematical model and developing a rapid heat sink evaluation method by making some mathematical treatments, that is, the implicit finite difference method (IFDM), to a heat conduction equation in the spherical coordinate system for the LED heat sink. With this method, we can estimate the static heat distribution and the temperature in different designs of LED heat sink. It can be more quickly achieved than the explicit finite difference method (EFDM) because it is unconditionally stable.The paper is organized as follows. In Section2, modelling of LED heat sink is introduced. In Section 3, we propose the discretization schemes and give some mathematical treatments. Some numerical simulations of heat distribution for the LED heat sink are provided in Section 4, and the conclusions are made in Section 5.


## 2. Modelling for LED Heat Sink

To the commonly used LED lamp (Figure1), regardless of its top cover and base, heat conduction equation for the heat sink can be expressed in the Cartesian coordinate system as follows:(1)∂U∂t=α·∂2U∂x2+∂2U∂y2+∂2U∂z2+qvρc·δx,y,z,where(2)δx,y,z=1,x,y,z=x0,y0,z00,else.We suppose that the LED lamp is continuously homogeneous and isotropic, while all variables are with the international standard unit system, whose domains are also supposed to be workable. Heat source lies near the top of LED lamp. For simplicity, we suppose that there is only one heat source within the LED lamp in this paper, which lies at (x0,y0,z0). Besides the top heat source, the inside longitudinal direction figure of LED heat sink can be described in Figure 2. The blank place is air, while in general the red is aluminum.Figure 1
LED lamp.Figure 2
LED inside longitudinal direction figure.Concerning the influence of heat source to the LED lamp, the popular treatment is to restrict the influence at slender grids, which closely surround(x0,y0,z0) in very limited grids. Obviously, it does not so well conform to the reality. The influence of heat source on LED should be continuous instead of being truncated by the limited grids. More reasonably, we can describe the influence with the pulse function; then the last part of (1) can be written as (3)qvρc·exp⁡-Ax-x02+y-y02+z-z02L2.L represents the height of LED lamp. A is a constant which is decided by the specific LED shape and power, which represents the influence of heat source to each part of LED lamp. We let a=A/L2 simplify (3), where the unit of a is m-2. Then the heat conduction equation (1) becomes(4)∂U∂t=α∂2U∂x2+∂2U∂y2+∂2U∂z2+qvρce-ax-x02+y-y02+z-z02.In our following work, with convenience, we consider the top and the base of heat sink as spherical instead of reality planar, and the position of heat source is limited at (0,0,z0) (z0=max⁡(z)). All segmentation is supposed equally, so the number of surrounding fins of the lamp body can be chosen as an even number. In considering spherical coordinate system, the heat conduction equation (4) becomes(5)∂U∂t=α∂2U∂r2+1r2∂2U∂θ2+1r2sin2θ∂2U∂φ2+2r∂U∂r+1r2tan⁡θ∂U∂θ+qvρce-ar2-2rr^cos⁡θ+r^2,where r^0≤r≤r^, 0≤θ<θ^1, 0≤φ≤2π, and 0≤t≤tf, or r^0≤r≤r^,θ^1≤θ≤θ^2, 2I-22π/N≤φ≤2I-12π/N, 0≤t≤tf, and I=1,2,…,N/2. (0,0,r^) is the location of heat source. N/2 represents the total number of surrounding fins of the lamp body (N is even). α, ρ, c, and κ are all piecewise constants, which are defined as follows.(I) Whenr=r^0 or r=r^, 0≤θ<θ^1, 0≤φ≤2π; θ^1≤θ≤θ^2, 2I-22π/N≤φ≤2I-12π/N, I=1,2,…,N/2, we have α=α2, κ=κ2, ρ=ρ2, c=c2.(II) Whenr^0<r<r^, 0≤θ≤θ^0(θ^0<θ^1), 0≤φ≤2π, we have α=α1, κ=κ1, ρ=ρ1, c=c1.(III) Whenr^0<r<r^, θ^0<θ<θ^1, 0≤φ≤2π; θ^1≤θ≤θ^2, 2I-22π/N≤φ≤2I-12π/N, I=1,2,…,N/2, we have α=α2, κ=κ2, ρ=ρ2, c=c2.We consider the relevant boundary conditions as follows.(IV) Whent=0, r^0≤r≤r^, 0≤θ<θ^1, 0≤φ≤2π or θ^1≤θ≤θ^2, 2I-22π/N≤φ≤2I-12π/N, I=1,2,…,N/2, we have(6)Ur,θ,φ=Uamb,where Uamb is a constant.(V) Whenr=r^0, 0≤θ<θ^1, 0≤φ≤2π, or θ^1≤θ≤θ^2, 2I-22π/N≤φ≤2I-12π/N, I=1,2,…,N/2, we have(7)-κ2∂U∂rr=r^0=H1Uamb-Ur=r^0.(VI) Whenr=r^, 0≤θ<θ^1, 0≤φ≤2π or θ^1≤θ≤θ^2, 2I-22π/N≤φ≤2I-12π/N, I=1,2,…,N/2, we have(8)κ2∂U∂rr=r^=H1Uamb-Ur=r^.(VII) Whenθ=0, r^0≤r≤r^, 0≤φ≤2π, we have(9)∂U∂θθ=0=0.(VIII) Whenθ=θ^0, r^0≤r≤r^, 0≤φ≤2π, we have(10)U1θ=θ^0=U2θ=θ^0,κ1∂U1∂n→θ=θ^0=κ2∂U2∂n→θ=θ^0,where U1=Uθ≤θ^0 and U2=Uθ≥θ^0.(IX) Whenθ=θ^1, r^0≤r≤r^, 2I-12π/N≤φ≤2I2π/N, I=1,2,…,N/2, we have(11)κ2∂U∂θθ=θ^1=H2Uamb-Uθ=θ^1.(X) Whenφ=2I-22π/N or φ=2I-12π/N, r^0≤r≤r^, θ^1≤θ≤θ^2, I=1,2,…,N/2, we have(12)-κ2∂U∂φφ=2I-22π/N=H2Uamb-Uφ=2I-22π/N,κ2∂U∂φφ=2I-12π/N=H2Uamb-Uφ=2I-12π/N.(XI) Whenθ=θ^2, r^0≤r≤r^, 2I-22π/N≤φ≤2I-12π/N, I=1,2,…,N/2, we have (13)κ2∂U∂θθ=θ^2=H2Uamb-Uθ=θ^2.Specifically, in the following work, we takeN=8 as an example. Then the longitudinal and crosscutting map of the LED heat sink considered in this paper could be shown in Figures 3 and 4. The color bar which lies in the right side of each map represents temperature. That is to say that we only consider four surrounding fins in our work. Other numerical values of N will be able to be finished similarly.Figure 3
The longitudinal map of the LED heat sink considered in this paper.Figure 4
The crosscutting map of the LED heat sink considered in this paper.


## 3. Discretization Schemes and Mathematical Treatments

Now we process the discretization at the points(tn,ri,θj,φm) for the heat equation (5), where(14)tn=nτ,τ=tfN1,n=0,1,…,N1;ri=r^0+i-0.5h1,h1=r^-r^0N2,i=1,2,…,N2;θj=j-0.5l1,l1=θ^0M1,j=1,2,…,M1;θj=θ^0+j-M1-0.5l2,l2=θ^2-θ^0M3-M1,j=M1+1,…,M2,…,M3;θ^1=θ^0+M2-M1l2.When0≤θ<θ^1 (i.e., j=1,2,…,M2), we let(15)φm=mh2,h2=2πN3,m=1,2,…,N3.Whenθ^1≤θ≤θ^2 (i.e., j=M2+1,…,M3), we let(16)φm=m-0.5h3,h3=π4N4,where m=1,2,…,N4,2N4+1,…,3N4,4N4+1,…,5N4,6N4+1,…,7N4.For simplicity, we denoteU(tn,ri,θj,φm) as Un,i,j,m afterwards. For the IFDM used in our paper, the partial derivative about t is approximated with the forward difference scheme, while the others are approximated with the central difference schemes. Particularly, we only indicate the treatments about φ and r.For the periodicity ofφ, when 0≤θ<θ^1 (i.e., j=1,2,…,M2), we let the related central differences be the following forms.Whenm=1, n=0,1,…,N1, i=1,2,…,N2, and j=1,2,…,M2, we let(17)∂2U∂φ2≈Un+1,i,j,2-2Un+1,i,j,1+Un+1,i,j,N3h22.Whenm=N3, n=0,1,…,N1, i=1,2,…,N2, and j=1,2,…,M2, we let(18)∂2U∂φ2≈Un+1,i,j,1-2Un+1,i,j,N3+Un+1,i,j,N3-1h22.For the central differences aboutr, we only show the following case.Wheni=1, n=0,1,…,N1 and j=1,2,…,M2, m=1,2,…,N3 or j=M2+1,…,M3, m=1,…,N4, 2N4+1,…,3N4,4N4+1,…,5N4,6N4+1,…,7N4, we let(19)∂U∂r≈Un+1,2,j,m-U^n+1,1,j,m2h1,∂2U∂r2≈Un+1,2,j,m-2Un+1,1,j,m+U^n+1,1,j,mh12,where U^n+1,1,j,m is the virtual symmetrical point of Un+1,1,j,m about r=r^0. With (7), we can get(20)-κ2Un+1,1,j,m-U^n+1,1,j,mh1=H1Uamb-Ur=r^0.We make the approximation(21)Ur=r^0≈U^n+1,1,j,m+Un+1,1,j,m2.So U^n+1,1,j,m can be deduced, and ∂U/∂r and ∂2U/∂r2 in this case are solved.Similarly, with (8)–(13), other partial derivatives can be finished. With these discretization schemes, the truncation error of (5) is O(τ+h12+l2+h^2), where l=max⁡(l1,l2) and h^=max⁡(h2,h3). With denser discretization, the result will turn better.Now, after the discretization of all variables, we can arrange the discretization of (5) as(22)Un=MUn+1+V,n=0,1,…,N1,where Un is (N2M2N3+4N2N4(M3-M2))×1 matrix, whose elements are Un,i,j,m, i=1,2,…,N2, and j=1,2,…,M2, m=1,2,…,N3, or j=M2+1,…,M3, m=1,2,…,N4,2N4+1,…,3N4,4N4+1,…,5N4,6N4+1,…,7N4. The point sequences (i,j,m) can be ordered as 1,j′,m′, 1,j′′,m′′, (2,1,1), …, (N2,M3,7N4), where j′=1,2,…,M2, m′=1,2,…,N3, j′′=M2+1,…,M3, and m′′=1,…,N4,2N4+1,…,3N4,4N4+1,…,5N4, 6N4+1,…,7N4. All of them are sorted in ascending order. U0 is (N2M2N3+4N2N4(M3-M2))×1 matrix, whose every element is Uamb. M is (N2M2N3+4N2N4(M3-M2))×(N2M2N3+4N2N4(M3-M2)) block matrix, whose every block is similar to a band matrix. But the band characteristic is destroyed by the disordered elements in the first row and last row. V is also (N2M2N3+4N2N4(M3-M2))×1 matrix, which greatly relies on qv.To solve the heat equation (5), we should run iteration (22) step by step. This treatment is unconditionally stable for the implicitness [24]. Our method is independent of the discretization step about heat transfer time (τ), so it can be rapidly achieved. It is very convenient and fast.


## 4. Numerical Simulations of Heat Distribution

For all heat distribution simulations here, we take a 0.5 w LED lamp as an example, where only a quarter of height is considered. For its shape structure, we letθ^0=53π/1000, θ^1=3π/50, θ^2=81π/1000, r^0=11/125m, r^=99/1000m, ρ1=1.205kg/m3, ρ2=2.7×103kg/m3, c1=1030J/(kg·K), c2=900J/(kg·K), κ1=0.026W/(m·K), and κ2=237W/(m·K). The ambient temperature Uamb is 25.00∘C. The numerical value of H1(H2) is fixed at 1.5. The radius R and height h of heat source are both 2 mm, so qv=1.99×107W/m3. To obtain a, we just consider e-a(r2-2rr^cos⁡θ+r^2). It should satisfy e-a(r2-2rr^cos⁡θ+r^2)→1 when (r,θ)→(r^,0), while R→2mm and h→2mm at the same time. After a short Matlab program, a=6.7×105m-2 approximately. The heat source influence is shown in Figure 5.Figure 5
Heat source influence.Numerical Simulation 1. To prove the validity of our presented schemes, we let qv=0W/m3 as verification. For the discretization of all variables, we let N2=14, M1=6, M2=8, M3=14, N3=16, and N4=8, and the heat transfer time is tf=3min. When N1=36, after a few seconds of Matlab program, the stable solution is represented that the final temperature of LED keeps 25.00∘C. When we prolong the heat transfer time, the result is still the same. It has no change since the heat transfer does not happen under this certain circumstance. This simulation indicates that our presented schemes are correct.Numerical Simulation 2. The discretization is the same as numerical simulation 1; that is, N2=14, M1=6, M2=8, M3=14, N3=16, N4=8, and the heat transfer time is tf=30 s. When N1=6 (τ=5s), the final temperature of LED ranges from 25.21∘C to 171.37∘C. This result is reasonable. That is the result for ignoring the heat dissipation through convection; in addition our built mathematical model is different from the actual model. We let Figure 6 describe the temperature distribution of the LED heat sink under this circumstance. In this paper, we only show figures of the temperature distributions at φ=0 for the symmetrical shape. When N1=15 (τ=2s), the final temperature of LED ranges from 25.19∘C to 171.36∘C. The temperature distribution is shown in Figure 7. From the two simulations, we can find that the calculation about time step τ has a little influence on the final result. Thus, in the next numerical simulations, we fix τ=5 s. In contrast, with the same discretization, the related explicit difference schemes need τ⩽0.001 s when we want to obtain the same result, so our present method is very fast and convenient.Figure 6
Temperature distribution attf=30 s and N1=6.Figure 7
Temperature distribution attf=30 s and N1=15.Numerical Simulation 3. When the discretization is N2=14, M1=6, M2=8, M3=14, N3=16, and N4=8, which is like the above two numerical simulations, we prolong the heat transfer time to tf=3min(N1=36); then we get the stable solution representing the final temperature ranges from 25.93∘C to 172.10∘C. The temperature distribution is shown in Figure 8. When the heat transfer time is tf=30min(N1=360), the stable solution represents the final temperature ranges from 26.88∘C to 173.06∘C. Its distribution of temperature is shown in Figure 9. From Figures 7, 8, and 9, we can find that the heat transfer becomes much clearer and better conformed to the reality when the heat transfer time tf becomes larger.Figure 8
Temperature distribution attf=3min.Figure 9
Temperature distribution attf=30min.Numerical Simulation 4. The discretization is denser than the above three numerical simulations, N2=18, M1=10, M2=12, M3=18, N3=16, and N4=8, when the heat transfer time is tf=3min(N1=36); then the stable solution represents the final temperature ranges from 26.10∘C to 275.49∘C, whose distribution of the temperature is shown in Figure 10, while the discretization becomes sparser than the above three numerical simulations, N2=8, M1=4, M2=6, M3=12, N3=8, and N4=4; the heat transfer time is still tf=3min⁡N1=36; then we get the stable solution representing the final temperature ranges from 25.12∘C to 34.77∘C. In this case, the related distribution of temperature is shown in Figure 11. Obviously, Figures 10 and 11 are far away from the reality for the too much dense or sparse mesh dissections.Figure 10
Temperature distribution attf=30 s (dense).Figure 11
Temperature distribution attf=30 s (sparse).


## 5. Conclusions

From the above numerical simulations of heat distribution, based on our fairly realistic mathematical model and treatments, we can design the unconditionally stable scheme of the heat distribution for LED heat sink. It is effective and convenient to evaluate the heat distribution in the spherical coordinate system. But the mesh dissections aboutr, θ, and φ should be chosen properly, because the too dense or sparse mesh dissections will become far away from the reality result. To the technique of computer, the large size of iterative matrix M which is a multiplication of mesh nodes about N2, M2, M3, N3, and N4 is a great challenge when we want much denser mesh dissections. In addition, when the size of M becomes larger, the evaluation speed becomes slower as a result. In spite of these shortcomings, our built mathematical model and unconditionally stable treatments still have advantage on the optimization of LED heat sink. With our treatments, the heat distribution and the temperature in different designs of LED lamp can be rapidly estimated. Furthermore, the built spherical coordinate system is better for the shape optimization design than in the Cartesian coordinate system, for the variable in the spherical coordinate system is in great accordance with actual angle shape. It is promising for the optimization design and actual engineering application of LED heat sink.



---

*Source: 3468246-2016-11-10.xml*