OPEN ACCESS
The aim of this study is to analyze a single component heat sink design numerically and analytically. The analysis included various modeling assumptions: lumped parameter; one dimensional and two dimensional models. It is shown that the lumped model for the considered problem is capable to accurately predict the steady state temperature of the heat sink with a small discrepancy at the center of the heat sink. Some results are presented graphically for various parameters to show their effect on the steady state temperature of the heat sink including power input; convective heat transfer coefficient; ambient temperature; thickness of the heat sink and the shape of the cross section of the heat sink (circular, ellipsoidal, squared). A few particular applications include single solar cell under high concentration; thermoelectric cooling device and/or highpower led. The findings of the research are twofold: educational and the derived formulae serve as a quick estimate tool to predict the effect of the various parameters on the heat sink temperature and thus aids to estimate the thermodynamic efficiency of the system.
heat sink, passive cooling, heat dissipation, single component, solar cell, thermoelectric cooling, led
Heat sinks [110] are devices used in high power systems to remove large amount of heat generated by its components [7]. A heat sink is used to increase the surface area which dissipates the heat faster and keeps the components under safe operating temperature ranges according to the manufacturer specifications [7]. The heat sink dissipates heat to its surroundings (ambient) passively or actively; with several working conditions and materials [8]. The heat is absorbed into the heat sink material via heat conduction and the heat is dissipated to the ambient either by convection or radiation [8]. Heat sinks are designed using thermal conductive materials—like copper and aluminum—and they work by dissipating heat through liquid cooling, natural convection, forced convection, or radiation [14]. A heat sink transfers the thermal energy generated by an electronic assembly or any heat generating component into a cooling medium. The heat is transferred from a higher temperature region to lower temperature region (fluid medium gas or liquid) by conduction, convection, radiation or by a combination of these heat transfer modes [14]. The performance of the heat sink is determined by many factors including the velocity of the coolant fluid, the thermal conductivity of the material, the thermal interface material, and the attachment method [14].
Specific applications are in order. Their performance and optimization were based on simulation. The use of highpower light emitting diodes (HPLED) for public illumination is an emerging subject, triggered by recent developments of different technologies including semiconductor materials [58], fluorescence techniques [9], driver electronics [10] or thermal control [11] among others including concentrated solar cells [1215].
In the next paragraphs, some applications and design parameters are reviewed briefly, just to give a few.
Heat sink configuration is very essential in heat sink design and optimization. As an example, as was demonstrated in [6], an analytical simulation model has been developed for predicting and optimizing the thermal performance of bidirectional fin heat sinks in a partially confined configuration.
As was pointed in [7], the increase in dissipated power per unit area of electronic components sets higher demands on the performance of the heat sink. The heat sink in this case, is a device used in computers to remove the large amount of heat generated by components, including integrated circuitsICs such as central processing unitsCPUs, chipsets and graphics cards, during their operation. The heat sink is used to increase the surface area which dissipates the heat faster and keeps the ICs under safe operating temperature and this is achieved by carful design conditions.
Another important factor in heat sink design is related to the material used and its composition is described in [8]. It is shown that heat sinks made of composite materials containing nonmetallic constituents, with a thermal conductivity as much as an order of magnitude less than typical metallic heat sinks, can provide an effective alternative where performance, cost and manufacturability are of importance.
For high power density as was highlighted in [9], thermal management is very essential especially in systems with all components including capacitors, inductors and semiconductor devices. These components used in power converters have maximum operating temperatures defined by the manufacturer. Increase in power density is the main factor which influences the thermal management to become so important and have to be treated properly.
Thermoelectric (TE) devices as was discussed in [10] can provide clean energy conversion and are environmentally friendly. The suggested application of a thermoelectric generator is to exploit the natural temperature differences between the air and the soil to generate small amounts of electrical energy. Since the conversion efficiency of even the best thermoelectric generators available is very low, the performance of the heat sinks providing the heat flow is critical. An experimental investigation was conducted to evaluate the performance of a thermoelectric (TE) module fitted to a conventional fin heat sink with a similarly sized heat source. [10]
In a different study, highpowered LED chips were considered [11]. The study presented a numerical analysis and experimental measurements of the temperature stabilization of highpower LED chips that they have obtained by employing an aluminum passive heat sink, designed to be used in a compact light bulb configuration. The study demonstrated that their system [11] kept the temperature of the LED chip wellbelow 70◦C yielding longterm operation of the device. The simulations have been performed for a lowcost device ready to install in public streetlights. The experimental measurements performed in different configurations show a nice agreement with the numerical calculations [11].
Solar energy applications were considered for concentrated photovoltaic solar cells. A thermal model for concentrator solar cells based on energy conservation principles was designed as was described by [12]. The solar cell was under 400X concentration with no cooling aid, the cell temperature would get up to about 1200 °C. Metal plates were used as heat sinks for cooling the system, which remarkably reduces the cell temperature. For a fixed concentration ratio, the cell temperature reduced as the heat sink area increased. In order to keep the cell at a constant temperature, the heat sink area needs to increase linearly as a function of the concentration ratio. GaInP/GaAs/Ge triplejunction solar cells were fabricated to verify the model. A cell temperature of 37 °C was measured when using a heat sink at 400X concentration [12].
In a different present study [13], an analysis of the benefits of passive cooling for high concentrator photovoltaic (HCPV) systems in terms of costs and kWh annual energy yields were provided. The performance of the heat sink has been related to the calculated energy yield of a standard triplejunction GaInP/GaAs/Ge HCPV cell in a system deployed at several suitable locations across the globe. Copper and Aluminum heat sinks have been considered and their merits have been compared. The finite element analysis software package COMSOL was employed to gain insights regarding a simple flat plate heat sink [13].
According to [14], the heat sink geometry is another important factor. The effect of heat sink geometry was addressed in [14]. The heat transfer performance of various commonly used fin geometries was compared. Realistic, manufacturable geometries were considered for minimizing thermal resistance at moderate laminar air velocities and pressure gradients. These consist of plate fins or pin fins, which can be round, elliptical, or square. The plate fins can be continuous (parallel plates) or staggered. The pin fins can be inline or staggered arrays [14].
Modeling and Experimental Evaluation of a single component Passive Heat Sinks for Miniature HighFlux Photovoltaic Concentrators were presented in reference [15]. The analysis was performed numerically and analytically. The obtained analytical results included Bessel functions which requires more than intuitive skills and greater effort to perform the design.
To summarize, the importance of heat sink is addressed in the previous paragraphs and the different design parameters and applications were briefly reviewed.
The currents study aims to reconsider the single component heat sink design for concentrated solar cell applications that was considered in [1113, 15] and LED application [14] taking into account several design parameters in order to derive a simple working formula which is capable to demonstrate these effects, easily and with acceptable accuracy. The methods of analysis used in the study are numerical and analytical in their nature. These methods were used to analyze different models which include different heat transfer modes (conduction, convection and radiation). The numerical analysis is used to check and compare the approximate analytical formulas derived in the following sections.
The rest of the article is arranged as follows: the problem statement is given in section 2, materials and methods are described in section 3, lumped parameter models are addressed in section 4, parametric study and numerical examples are given in section 5, summary and conclusions are given in section 6, the numerical scheme of the complete model is given in Appendix I and finally the Microsoft excel VBA code is given in Appendix II.
Consider a heat sink to passively dissipate an amount of power p (w) absorbed at its center within an area within a diameter of size d0 (m). The heat sink cross section could be a circle (Figure 1a); an ellipse (Figure 1b) or a square (Figure 1c). The diameter of the heat sink is d (m) for the circular cross section.
Figure 1. Schematics of the heat sink a  circular cross section; b – ellipsoidal cross section; c – squared cross section
The dimensions of the ellipse and the square are determined such that the volume v (m^{3}) of the heat sink is fixed. The heat sink is in the shape of a plate with thickness z_{max} (m). The heat sink material is a metal with high thermal conductivity k (w/m/K) such as aluminum; copper; or iron. The heat energy is dissipated by heat convection and thermal radiation. The heat convection is modeled by using Newtonian heat transfer law with convective heat transfer coefficient h (w/m^{2}/K) and the thermal radiation is modeled by using StefanBoltzmann law. Heat from the heat sink is dissipated to the surrounding with fixed ambient temperature Ta K. The temperature distribution T K of the heat sink could be determined by using the energy equation and the appropriate boundary conditions. The rate of energy E J balance with respect to time t s is given by (all units are expressed using SI international system of units):
$\frac{\partial E}{\partial t}=\nabla \cdot(k \nabla T)$ (1)
After applying equation (1) for cylindrical coordinates with a symmetric shape (circular disk) and for constant thermal conductivity coefficient, the energy balance equation is writing in the following explicit form:
$\frac{\partial(\rho c T)}{\partial t}=k\left(\frac{\partial^{2} T}{\partial r^{2}}+\frac{1}{r} \frac{\partial T}{\partial r}\right)+k \frac{\partial^{2} T}{\partial z^{2}}$ (2)
where $\rho$ kg/m^{3} is the density and c J/kg/K is the heat capacity of the heat sinks' material (solid) and r m is the radius of the heat sink and z m is the vertical distance along the height of the disk (cylinder).
After some mathematical manipulation and assuming constant thermophysical properties of the materials, the energy balance equation is written as follows [1617]:
$\frac{1}{\alpha} \frac{\partial T}{\partial t}=\left(\frac{\partial^{2} T}{\partial r^{2}}+\frac{1}{r} \frac{\partial T}{\partial r}\right)+\frac{\partial^{2} T}{\partial z^{2}}$ (3)
where $\alpha$ is the thermal diffusivity ($\alpha=\frac{k}{\rho c}$).
At zero radius Eq. (3) is replaced by $\frac{1}{\alpha} \frac{\partial T}{\partial t}=2 \frac{\partial^{2} T}{\partial r^{2}}+\frac{\partial^{2} T}{\partial z^{2}}$.
This equation is linear one dimensional partial differential equation in time and two dimensional in space. In order to fully specify the aforementioned temperature distribution, initial temperature condition is required Ti along with four boundary conditions at the edges of the heat sink. The boundary conditions are given as follows:
At r=0,
$\frac{\partial T}{\partial r}=0$ (4)
At r=R (=d/2 radius of the heat sink),
$k \frac{\partial T}{\partial r}=h\left(TT_{a}\right)+\sigma\left(T^{4}T^{4}\right)$ (5)
At z=0,
$\left\{\begin{array}{ll}{k \frac{\partial T}{\partial z}+h\left(TT_{a}\right)+\sigma\left(T^{4}T^{4}\right)=p} & {r \leq \frac{d_{0}}{2}} \\ {k \frac{\partial T}{\partial z}+h\left(TT_{a}\right)+\sigma\left(T^{4}T^{4}\right)=0} & {r>\frac{d_{0}}{2}}\end{array}\right.$ (6)
And at z=z_{max},
$k \frac{\partial T}{\partial z}=h\left(TT_{a}\right)+\sigma\left(T^{4}T_{a}^{4}\right)$ (7)
Eq. (3) alongside with the boundary conditions (Eqns. 47) constructs the mathematical model of the heat sink. This set of equations is nonlinear due to the nonlinearity of some of the boundary conditions, its solutions could be achieved numerically (the numerical method of solution will be discussed in the next section). Fortunately, the high conductivity of the solid material will simplify the partial differential equation from two dimensional to a lumped model.
As mentioned earlier at the end of the previous section, the solution of the mathematical model could be achieved numerically (details of the numerical solution will be presented in appendix I – the numerical scheme and in appendix II – its implementation using Microsoft Excel VBA). The materials of the heat sink analyzed in this study include aluminum, copper and iron. The following table (table 1) summarizes the properties used in the numerical examples section.
Table 1. Thermophysical properties of the heat sink materials used in the study
Material 
Thermal conductivity w/m/K 
Density Kg/m^{3} 
Heat capacity J/kg/K 
Thermal diffusivity m^{2}/s 
Aluminum type 6061T6 
167 
2700 
902 
6.8572x10^{5} 
Copper 
401 
8960 
385 
11.6250x10^{5} 
Iron 
80.4 
7874 
449 
2.2741x10^{5} 
The geometry dimensions (diameter and thickness) are picked such that the volume is fixed to 25cc [15]. Table 2 summarizes the value of the parameters and constant used to produce the plots given in the examples section.
The combination of the high thermal conductivity of the solid materials and the relatively small characteristic length L_{c} with passive cooling, facilitate the simplification of the mathematical model by means of the Biot number ($=\frac{h L_{c}}{k}$) [18]. The temperature distribution of the heat sink could be considered uniform for Biot numbers less than 0.1 [19]. For the vertical direction Biot_{z} =0.00015 and for the radial direction Biot_{r}=0.0012 (see table 2). These values suggest that the temperature of the heat sink could be considered uniform (lumped). This will be treated in the next section.
Table 2. Geometry dimensions and other parameters and constants used in the calculations
Name 
Value 
Diameter of the heat sink m^{3} 
0.08 
Thickness of the heat sink  m 
0.005 
Surrounding (ambient) temperature  K 
300 
Convective heat transfer coefficient w/m^{2}/K 
5 
Biot number =hz_{max}/k vertical Biot number =hd/2/k radial 
0.00015 0.0012 
Stefan – Boltzmann constant  w/m^{4}/K 

p=4*Arc tan(1) 
3.141592654 
In this section, the mathematical model will be simplified following the lumped system approach. Four cases will be considered: the complete lumped model and three approximations: heat convection dominated; thermal radiation dominated and approximate combined convectionradiation model.
4.1 Lumped heat sink model
Consider a heat sink as was introduced in section 2 (see Fig, 1) but now the temperature is uniform all over the space. The power absorbed by the heat sink is dissipated by modes of heat convection and thermal radiation heat transfer laws. The energy balance equation is given by:
$\frac{d(\rho c v T)}{d t}=ph A\left(TT_{a}\right)\sigma A\left(T^{4}T_{a}^{4}\right)$ (8)
And the heat transfer area A is given by:
$A=\frac{\pi d^{2}}{2}+\pi d z_{\max }$ (9)
Eq. (9) is nonlinear with respect to temperature and could be solved by means of RungeKutta method of order fourth (RK4). In the following subsections, further simplifications allow approximate analytical solutions to predict the temperature of the heat sink as a function of time.
4.2 Approximate convectionradiation model
Eq. (8) could be simplified by introducing the following substitution:
$T(t)=(1+\theta(t)) T_{a}$ (10)
where $\theta(t)$ is the normalized shift of the heat sink temperature relative to the ambient temperature. By introducing Eq. (10) into Eq. (8) and after mathematical manipulations, the following first order ordinary differential equation for the normalized heat temperature shift of the heat sink is given by:
$\tau \frac{d \theta}{d t}=\theta^{2}+c_{1} \thetac_{2}$ (11)
where $\tau=\frac{\rho c v}{6 \sigma A T_{a}^{3}}$; $c_{1}=\frac{2}{3}+\frac{h}{6 \sigma T_{a}^{3}}$; and $c_{2}=\frac{p}{6 \sigma T_{a}^{3}}$.
The solution of Eq. (11) is given by the following parametric relations:
$u(t)=c_{3} \frac{1+c_{4} e^{2 c_{1} \frac{t}{\tau}}}{1c_{4} e^{2 c_{3} \frac{t}{\tau}}} ; u(t)=\theta(t)+\frac{c_{1}}{2} ; c_{3}=\sqrt{c_{2}+\frac{c_{1}^{2}}{4}} ; c_{4}=\frac{u_{0}c_{3}}{u_{0}+c_{3}}$ (12)
The steady state heat sink temperature $\theta_{s}$ could be evaluated by the following expression:
$\theta_{s}=\frac{2^{c_{2}} / c_{1}}{\sqrt{1+4^{c_{2}} / c_{1}}+1}$ (13)
4.3 Heat convection model
Eq. (8) could be reduced to linear form after introducing the radiative heat transfer coefficient $h_{r}=\sigma\left(T^{2}+T_{a}^{2}\right)\left(T+T_{a}\right)$. The average temperature is used to estimate the value of h_{r}. The linear form of Eq. (8) is given by:
$\rho c v \frac{d T}{d t}=pA\left(h+h_{r}\right)\left(TT_{a}\right)=pA h\left(1+\xi_{c}\right)\left(TT_{a}\right)$ (14)
where $\xi_{c}$ is the ratio between the radiative heat transfer coefficient and the convective heat transfer coefficient. The solution of Eq. (14) is given by:
$T(t)=T_{a}+\frac{p}{A h\left(1+\xi_{c}\right)}\left(1e^{\frac{A h\left(1+\xi_{c}\right)}{\rho c v} t}\right)$ (15)
4.4 Heat radiation model
Eq. (8) is rewritten after introducing the ration between the convective heat transfer coefficient and the radiative heat transfer coefficient $\xi_{r}$ as follows:
$\rho c v \frac{d T}{d t}=pA\left(1+\xi_{r}\right) \sigma\left(T^{4}T_{a}^{4}\right)$ (16)
For convenience of the reader, Eq. (16) is rearranged after separation of variables as follows:
$\frac{d T}{\beta^{4}T^{4}}=\frac{d t}{2 \beta^{3} \tau_{r}}$ (17)
where $\beta=\sqrt{T_{a}^{4}+\frac{p}{A \sigma\left(1+\xi_{r}\right)}}$ is the steady state temperature and $\tau_{r}=\frac{\rho c v}{2 \beta^{3} A \sigma(1+\xi r)}$ is a time constant of the problem. After using the partial fractions decomposition of Eq. (17) followed by integrating the appropriate fractions, the final solution of Eq. (17) is given by:
$\frac{1}{2} \ln \left(\left(\frac{\beta+T}{\betaT}\right)\left(\frac{\betaT_{a}}{\beta+T_{a}}\right)\right)+\operatorname{Arctan}\left(\frac{T}{\beta}\right)\operatorname{Arctan}\left(\frac{T_{a}}{\beta}\right)=\frac{t}{\tau_{r}}$ (18)
Eq. (18) is an implicit relation of the heat sink temperature with respect to time.
5.1 Parametric study
The heat sink's basic configuration is a circular disk with a diameter of size 0.08 m and thickness of 0.005 m. the heat sink's material is aluminum. The volume of the heat sink is approximately 25 cc. The steady state temperature of the heat sink is given by Eq. (13) which is used to study its changes with respect to the following parameters:
Distributed vs. limped models: the two dimensional model is solved numerically followed by lumped models. The lumped models are motivated by small Biot numbers in the radial and axial directions of the disk heat sink.
Materials of the heat sink: the heat sink is made of metals with high conductivity coefficient. three materials are considered – aluminum, copper and iron.
Power input: the power input application dependentconcentrated solar cell, thermoelectric generator, or high power led. Values of the power input in the range 110 watts are used for demonstration.
Convective heat transfer coefficient: convective heat transfer ranges from natural convection to forced convection with different values of the convective heat transfer coefficient. Values in the range 110 w/m^{2}/K are used for demonstration.
Ambient temperature: the ambient temperature could be changed seasonally or geographically or by means of air conditioning according to specifications and requirements. Values of the ambient temperature in the range zero – 35 degrees Celsius are used for demonstration.
Thickness of the heat sink: the thickness of the considered heat sink attains approximately 10 % of the heat transfer to the ambient. Values in the range 15 mm are used for demonstration.
Geometric shape of the cross section of the heat sink: three options for the shape are consideredcircular, ellipsoidal and squared. The area for heat dissipation from the heat sink for the circular shape is given by Eq. (9). Similar expressions are given below for the ellipsoidal and squared shapes, under fixed volume condition. The volumes of the disk for circular shape; ellipsoidal shape; and the squared shape are given by:
$v=\pi \frac{d^{2}}{4} z_{\max } ; v=\pi a b z_{\max } ; \quad v=s^{2} z_{\max }$ (19)
where a, b are the major and minor radiuses of the ellipse and s is the side of the square.
For ellipsoidal shape Eq. (9) is replaced by:
$A=\pi \frac{d^{2}}{2}+\pi(a+b)\left(1+\frac{3 \delta}{10+\sqrt{43 \delta}}\right) z_{\mathrm{max}} ; \delta=\left(\frac{ab}{a+b}\right)^{2}$ (20)
In Eq. (20) the perimeter of the ellipse is calculated by Ramanujan approximation for circumference of an ellipse (https://www.johndcook.com/blog/2013/05/05/ramanujancircumferenceellipse/). In order to fully specify the dimensions of the ellipse, the minor radius b is arbitrarily assumed to be the ratio between the major parameter a and the golden section ratio $\phi=\frac{1+\sqrt{5}}{2}$, i.e. $\left(b=\frac{a}{\phi}\right)$.
Finally, for squared cross section, Eq. (9) is replaced by:
$A=\pi \frac{d^{2}}{2}+4 s z_{\max }$ (21)
5.2 Numerical examples
In this section, the analysis of the heat sink is demonstrated numerically and the results are presented graphically. The calculation is performed numerically and the partial differential equation (pde) is approximated with a forward finite difference scheme [16] and solved using Microsoft Excel's VBA. The calculation is stopped whenever steady state condition was reached, such that the dissipated heat to the ambient via convection and radiation equals the power input. The process of the solution was repeated for different materials (see Figure 2).
Figure 2. Dissipated power of the heat sink vs. time for different materials – aluminum (upper curve), copper (middle curve) and iron (lower curve). The input power was 2w
The steady state temperature distribution of the heat sink was observed for the considered materials as a function of the radial position measured from the center of the disk. The maximum temperature was observed at the center and it was slightly reduced towards the extreme of the disk. The temperature distribution is a function of the thermal conductivity with the following observation: the lower the thermal conductivity, the temperature (see Figure 3).
Figure 3. Steady state temperature distribution of the heat sink as a function of the radial position measured from the center of the disk. Three materials were considered: iron (upper curve), aluminum (middle curve) and copper (lower curve)
The temperature distribution suggested a uniform temperature of the heat sink and the distributed model was replaced by a lumped parameter model (Eq. 8). This equation was solved numerically by means of RK4 method and was approximated analytically (Eq. 12) and the results were plotted as a function of time (see Figure 4). Excellent agreement between the two results was observed.
The lumped parameter model was further simplified using convective model (Eq. 15) and radiative model (Eq. 18). The results are presented graphically (see Figure 5).
Figure 4. The temperature of the heat sink vs. time as was produced numerically by means of RK4 method which solved Eq. 8 and approximated analytically which is given by Eq. 12
Figure 5. Temperature of the heat sink of the convective model and the radiative model vs. time
Figure 6. Steady state temperature as a function of power input
The steady state temperature of the heat sink (Eq. 13) is studied parametrically to check for the effects of: power input (Figure 6); heat convection coefficient (Figure 7); ambient temperature (Figure 8); and thickness of the heat sink for different materials (Figure 9).
Effect of the power input on the steady state temperature of the heat sink is shown graphically (see Figure 6).
Effect of the convective heat transfer coefficient on the steady state temperature of the heat sink is shown graphically (see Figure 7).
Figure 7. Steady state temperature as a function of convective heat transfer coefficient
Effect of the ambient temperature on the steady state temperature of the heat sink is shown graphically (see Figure 8).
Figure 8. Steady state temperature as a function of ambient temperature
Figure 9. Steady state temperature as a function of the thickness of the heat sink. Different shapes of the heat sink were considered
Effect of the thickness of the heat sink on the steady state temperature of the heat sink is shown graphically (see Figure 9).
In this study, a single component heat sink was considered. The type of applications includes high concentration solar cells, thermoelectric power generation, high power LED and more.
The two dimensional heat sink model were solved numerically and its solution were compared with the lumped parameter model. The two dimensional model allows parametric studies and the effect of different materials (aluminum, copper and iron) was obtained.
The lumped parameter model was solved using the RK4 method and analytically using different approximations and modeling assumptions: convective, radiative and mixed. It was shown that the steady state temperature was in a good agreement while comparing the results achieved from the solution of the two dimensional model and from the lumped parameter model. This agreement is in line with the check of Biot numbers. The values of Biot numbers calculated radially and axially were less than 0.1 which means that lumped parameter assumption should be followed. In fact, small discrepancies exist especially near the center of the heat sink.
The good agreement of the steady state temperature achieved from the two dimensional model and lumped parameter model allowed the derivation of simple formulae to predict the steady state temperature of the heat sink and allowed parametric study of the several parameters involved in the problem. These parameters included: power input; convective heat transfer coefficient; ambient temperature; the thickness if the heat sink; and the effect of the geometrical shape of the heat sink (circular, ellipsoidal and squared). These effects were shown graphically (see Figures, 69).
The transient response was also achieved and the response time of the considered problem was approximately 30 minutes.
The proposed formulae for the steady state temperature of the heat sink could be used as a quick tool of the heat sink optimization and parametric calculations of the several involved parameters could be performed easily.
Finally, we recommend extending the results of this study by considering the application of single solar cell and its heat sink design for space applications, were air currents are missing, both experimentally and analytically.
A 
area, m^{2} 
C 
specific heat, J. kg^{1}. K^{1} 
E R T T_{a} c_{1}, c_{2} d d_{0} h k p 
energy, J radius of heat sink, m temperature of heat sink, ambient temperature, parameters used in equations 11, 12 diameter of heat sink, m diameter of irradiated area, m coefficient of heat convection, W m^{2} K^{1} thermal conductivity, W.m^{1}. K^{1} power W 
R z t 
radial direction, m axial direction, m time, s 
Greek symbols 

$\alpha$ 
thermal diffusivity, m^{2}. s^{1} 
$\beta$ 
parameter used in equation 17, 18 
$\rho$ 
desity, kg m^{3} 
Ɵ Ɵ_{s} 
dimensionless temperature used i equation 10 dimensionless temperature used i equation 13 
$\sigma$ $\tau$ 
Stefan Boltzmann coefficient, W. m^{2}. K^{4} time constant. S 
r * 
Radiation dimensionless temperature 
Subscripts n Index ir, iz 
time step
radial step, axial step 
APPENDIX I
Finite difference scheme of the heat sink model
The first order differential of the temperature with respect to time is approximated with forward difference scheme, but for radial and axial positions, a central difference schemes were used for the first and second differentials. The boundary conditions were treated first and the pde were treated second, to reflect the order of the numerical calculations while implemented using Microsoft Excel VBA. The radial position is divided into nr divisions and the thickness of the heat sink is divided into nz divisions. The radial position is defined by using ir index and the axial position is defined by using iz index. Both indices start at zero. The following equations are the finite difference version of the mathematical model.
Boundary condition at zero radius and zero height (increasing height faces downwards) is given by:
$p=k \frac{\left(T^{n}(0,1)T^{n}(0,1)\right)}{2 d z}+h\left(T^{n}(0,0)T_{a}\right)+\sigma\left(T^{n4}(0,0)T_{a}^{4}\right)$ (I.1)
where the superscript n denotes time step and the (i,j) indices denote radial and axial positions respectively. Eq, (I.1) is used for zero height position and all radiuses except zero.
Similarly, the boundary condition at maximum height is given by:
$k \frac{\left(T^{n}(i r, n z+1)T^{n}(i r, n z1)\right)}{2 d z} = h\left(T^{n}(i r, n z)T_{a}\right)+\sigma\left(T^{n 4}(i r, n z)T_{a}^{4}\right)$ (I.2)
Following the same procedure, the boundary condition at zero radius is given by:
$T^{n}(1, i z)=T^{n}(+1, i z)$ (I.3)
Similarly, the boundary condition at maximal radius is given by:
$k \frac{\left(T^{n}(n r+1, i z)T^{n}(n r+1, i z)\right)}{2 d r} = h\left(T^{n}(n r, i z)T_{a}\right)+\sigma\left(T^{n 4}(n r ; i z)T_{a}^{4}\right)$ (I.4)
The pde at zero radius is given by:
$\frac{1}{\alpha} \frac{\left(T^{n+1}(i r, i z)T^{n}(i r, i z)\right)}{\Delta t} = 2 \frac{T^{n}(i r+1, i z)2 T^{n}(i r, i z)+T^{n}(i r1, i z)}{\Delta r^{2}}+\frac{T^{n}(i r, i z+1)2 T^{n}(i r, i z)+T^{n}(i r, i z1)}{\Delta z^{2}}$ (I.5)
And finally, the pde everywhere except at zero radius is given by:
$\frac{1}{\alpha} \frac{\left(T^{n+1}(i r, i z)T^{n}(i r, i z)\right)}{\Delta t}=\frac{T^{n}(i r+1, i z)2 T^{n}(i r, i z)+T^{n}(i r1, i z)}{\Delta r^{2}}+\frac{\left(T^{n}(i r+1, i z)T^{n}(i r1, i z)\right)}{i r \Delta r 2 \Delta r}+\frac{T^{n}(i r, i z+1)2 T^{n}(i r, i z)+T^{n}(i r, i z1)}{\Delta z^{2}}$ (I.6)
Convergence and stability are granted for proper choice of time step [20].
APPENDIX II
Microsoft Excel VBA code
Sub sol()
Dim tre(1 To 9, 1 To 3) As Double
'temperature array 8x2 radius x thickness
'pproximately uniform in thickness
Pi = 4# * Atn(1#)
Cells(1, 4) = Pi
h = 5#
sig = 0.00000005667
d = 0.08
d0 = 0.005
Z = 0.005
dr = d / 2# / 8#
dz = Z / 2#
'aluminum type 6061T6
ro = 2700#
c = 902#
k = 167#
'
'copper
'ro = 8960#
'c = 385#
'k = 401#
'iron
'r = 7874#
'c = 449#
'k = 80.4
alfa = k / ro / c
dt = 0.5 / alfa / (1# / dr ^ 2 + 1# / dz ^ 2)
'resultcheck dt=0.036, take dt/10 for accuracy.
'pick dt=0.004 (250 time steps = 1 sec)
dt = 0.004
p = 2#
pdis = p * 4# / d0 ^ 2 / Pi
nz = 2
nr = 8
dtr = alfa * dt / dr ^ 2
dtz = alfa * dt / dz ^ 2
Ta = 300#
'initial state
For iz = 1 To nz + 1
For ir = 1 To nr + 1
tre(ir, iz) = Ta
Next ir
Next iz
'reset heat transfer
q1 = 0#
q2 = 0#
q3 = 0#
q = 0#
Timee = 0#
'loop over time
For it = 0 To 1000000
Timee = it * dt
'print state every second
If ((it Mod 250) = 0) Then
Cells(2 + it / 250, 14) = Timee / 60#
Cells(2 + it / 250, 15) = q
Cells(1, 3) = it
For iz = 1 To nz + 1
For ir = 1 To nr + 1
Cells(3 + iz, 4 + ir) = tre(ir, iz)
Next ir
Next iz
End If
'reset heat transfer
q1 = 0#
q2 = 0#
q3 = 0#
q = 0#
'r=0, z=0
tre(0, 1) = tre(0, 1) + 2# * dz / k * (pdis  h * (tre(0, 0)  Ta)  sig * (tre(0, 0) ^ 4  Ta ^ 4))
'r=0, z=0.005
tre(0, nz + 1) = tre(0, nz  1) + 2# * dz / k * (0# * pdis  h * (tre(0, nz)  Ta)  sig * (tre(0, nz) ^ 4  Ta ^ 4))
'bc at r=0, all z
For iz = 0 To nz
tre(1, iz) = tre(1, iz)
Next iz
'bc at r=R, all z
For iz = 0 To nz
tre(nr + 1, iz) = tre(nr  1, iz) + 2# * dr / k * (0# * pdis  h * (tre(nr, iz)  Ta)  sig * (tre(nr, iz) ^ 4  Ta ^ 4))
Next iz
'z=0, all r
For ir = 1 To nr
tre(ir, 1) = tre(ir, 1) + 2# * dz / k * (0# * pdis  h * (tre(ir, 0)  Ta)  sig * (tre(ir, 0) ^ 4  Ta ^ 4))
Next ir
'z=h, all r
For ir = 1 To nr
tre(ir, nz + 1) = tre(ir, nz  1) + 2# * dz / k * (0# * pdis  h * (tre(ir, nz)  Ta)  sig * (tre(ir, nz) ^ 4  Ta ^ 4))
Next ir
'pde at r=0, all z
For iz = 0 To nz
tre(0, iz) = tre(0, iz) * (1#  4# * dtr  2# * dtz) + 4# * tre(1, iz) * dtr + dtz * (tre(0, iz + 1) + tre(0, iz  1))
Next iz
'pde elsewhere
For iz = 0 To nz
For ir = 1 To nr
tre(ir, iz) = tre(ir, iz) * (1#  2# * dtr  2# * dtz) + dtr * (tre(ir + 1, iz) + tre(ir  1, iz)) + dtr / 2# / ir * (tre(ir + 1, iz)  tre(ir  1, iz)) + dtz * (tre(ir, iz + 1) + tre(ir, iz  1))
Next ir
Next iz
For ir = 0 To nr  1
tav0 = (tre(ir, 0) + tre(ir + 1, 0)) / 2#
q1 = q1 + (h * (tav0  Ta) + sig * (tav0 ^ 4  Ta ^ 4)) * 2# * Pi * (ir + ir + 1) * dr / 2# * dr
tav1 = (tre(ir, nz) + tre(ir + 1, nz)) / 2#
q2 = q2 + (h * (tav0  Ta) + sig * (tav0 ^ 4  Ta ^ 4)) * 2# * Pi * (ir + ir + 1) * dr / 2# * dr
Next ir
For iz = 0 To nz  1
tavz = (tre(nr, iz) + tre(nr, iz + 1)) / 2#
q3 = q3 + (h * (tavz  Ta) + sig * (tavz ^ 4  Ta ^ 4)) * 2# * Pi * nr * dr * dz
Next iz
q = q1 + q2 + q3
Next it
End Sub.
[1] Kraus, A.D., BarCohen, A. (1995). Design and analysis of heat sinks, John Wiley and Sons, Inc., New York, (2nd printing 1997).
[2] BarCohen, A., Kraus, A.D. (1999). Wiley series in thermal management of microelectronics and electronic systems: Incorpera, F. P., Liquid cooling of electronic devices and singlephase convection. John Wiley and Sons, Inc., New York.
[3] BarCohen, A., Iyengar, M. (2002). Design and optimization of air cooled heat sinks for sustainable development. IEEE CPT Transactions, 25(4): 584591. https://doi.org/10.1109/TCAPT.2003.809112
[4] Iyengar, M., BarCohen, A. (2003). Least energy optimization of air cooled heat sinks for sustainable development. IEEE CPT Transactions, 25(4): 1625. https://doi.org/10.1109/tcapt.2003.811463
[5] Wang, C.C., Hung, C., Chen, W.H. (2012). Design of heat sink for improving the performance of thermoelectric generator using twostage optimization. Energy, 39: 236245. https://doi.org/10.1016/j.energy.2012.01.025
[6] Lee, S. (1995). Optimum design and selection of heat sinks. Eleventh IEEE SEMITHERMTM Symposium, 4854. https://doi.org/10.1109/95.477468
[7] Aravindh, B.S., Nair, T.R.G. (2008). Heat sink performance analysis through numerical technique. IEEE Symposium (NSSP08), Bangalore, pp. 17.
[8] Culham, J.R., Khan, W.A., Yovanovich, M.M., Muzychka, Y.S. (2001). The influence of material properties and spreading resistance in the thermal design of plate fin heat sinks. Journal of Electronic Packaging, 129(1): 7681. https://doi.org/10.1115/1.2429713
[9] Staliulionis, Z., Zhang, Z., Pittini, R., Andersen, M.A.E., Tarvydas, P., Noreika, A. (2014). Investigation of heat sink efficiency for electronic component cooling applications. Elektronika Ir Elektrotechnika, 20(1): 4954. https://doi.org/10.5755/j01.eee.20.1.6167
[10] Ong, K.S., Tan, C.F., Lai, K.C. (2017). Methodological considerations of using hermoelectrics with fin heat sinks for cooling applications. Applied Sciences, 7(2): 6272. https://doi.org/10.3390/app7020062
[11] Balvis, E., Bendaña, R., Michinel Alvarez, H., Fernández De Córdoba Castellá, PJ., Paredes, A. (2015). Analysis of a passive heat sink for temperature stabilization of highpower LED bulbs. Journal of Physics: Conference Series, 605: 01200510120058. https://doi.org/10.1088/17426596/605/1/012005
[12] Cui, M., Chen, N.F., Yang, X.L., Wang, Yu., Bai, Y.M., Zhang, X.W. (2009). Thermal analysis and test for single concentrator solar cells. Journal of Semiconductors, 30(4): 6164. https://doi.org/10.1088/1674.4926/30/4/044011
[13] Micheli, L., Fernández, E.F., Almonacid, F., Mallick, T.K., Smesta, G.P., (2016). Performance, limits and economic perspectives for passive cooling of High Concentrator Photovoltaics. Solar Energy Materials and Solar Cells, 153: 164178. https://doi.org/10.1016/j.solmat.2016.04.016
[14] Soodphakdee, D., Behnia, M., Copeland, D.W. (2001). A comparison of fin geometries for heat sinks in laminar forced Ccnvection: Part I  round, elliptical, and plate fins in staggered and inline configurations. The International Journal of Microcircuits and Electronic Packaging, 24(1): 6876. https://doi.org/10.1088/16744926/30/4/044011
[15] Sun, J., Israeli, T., Reddy, T., Scoles, K., Gordon, J.M., Feuermann, D. (2004). Modeling and experimental evaluation of passive heat sinks for miniature highflux photovoltaic concentrators. ASME International Solar Energy Conference, 431440. https://doi.org/10.1115/ISEC200465184.
[16] Smith, G.D. (1985). Numerical Solutions of Partial Differential Equations: Finite Difference Method. 3rd edition, Oxford University Press, NY, 75.
[17] Smith, G.D. (1985). Numerical Solutions of Partial Differential Equations: Finite Difference Method. 3rd edition, Oxford University Press, NY, 76.
[18] Incropera, F.P., DeWitt, D.P., Bergman, T.L., Lavine, A.S., (2007). Fundamentals of Heat and Mass Transfer. 6th edition, New York: J. Wiley, 260.
[19] Incropera, F.P., DeWitt, D.P., Bergman, T.L., Lavine, A.S., (2007). Fundamentals of Heat and Mass Transfer. 6th edition, New York: J. Wiley, 261.
[20] Smith, G.D. (1985). Numerical Solutions of Partial Differential Equations: Finite Difference Method。 3rd edition, Oxford University Press, NY, 43.