# Parametric Study of Uniform Heat Sink Design

Parametric Study of Uniform Heat Sink Design

Mahmoud M. Huleihil

The Arab Academic Institute of Education, Beit-Berl College, Kfar Saba 44905, Israel

Corresponding Author Email:
Mahmud.ana@gmail.com
Page:
201-209
|
DOI:
https://doi.org/10.18280/i2m.180216
12 January 2019
|
Accepted:
23 March 2019
|
Published:
30 June 2019
| Citation

OPEN ACCESS

Abstract:

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 high-power 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.

Keywords:

heat sink, passive cooling, heat dissipation, single component, solar cell, thermoelectric cooling, led

1. Introduction

Heat sinks [1-10] are devices used in high power systems to remove large amount of heat generated by its components . 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 . The heat sink dissipates heat to its surroundings (ambient) passively or actively; with several working conditions and materials . 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 . 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 [1-4]. 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 [1-4]. 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 [1-4].

Specific applications are in order. Their performance and optimization were based on simulation. The use of high-power light emitting diodes (HP-LED) for public illumination is an emerging subject, triggered by recent developments of different technologies including semiconductor materials [5-8], fluorescence techniques , driver electronics  or thermal control  among others including concentrated solar cells [12-15].

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 , 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 , 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 circuits-ICs such as central processing units-CPUs, 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 . It is shown that heat sinks made of composite materials containing non-metallic 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 , 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  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. 

In a different study, high-powered LED chips were considered .  The study presented a numerical analysis and experimental measurements of the temperature stabilization of high-power 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  kept the temperature of the LED chip well-below 70◦C yielding long-term operation of the device. The simulations have been performed for a low-cost device ready to install in public streetlights. The experimental measurements performed in different configurations show a nice agreement with the numerical calculations .

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 . 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 triple-junction solar cells were fabricated to verify the model. A cell temperature of 37 °C was measured when using a heat sink at 400X concentration .

In a different present study , 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 triple-junction 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 .

According to , the heat sink geometry is another important factor. The effect of heat sink geometry was addressed in . 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 .

Modeling and Experimental Evaluation of a single component Passive Heat Sinks for Miniature High-Flux Photovoltaic Concentrators were presented in reference . 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 [11-13, 15] and LED application  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.

2. Problem Statement

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.

## 1.png 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 (m3) of the heat sink is fixed. The heat sink is in the shape of a plate with thickness zmax (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/m2/K) and the thermal radiation is modeled by using Stefan-Boltzmann 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/m3 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 thermo-physical properties of the materials, the energy balance equation is written as follows [16-17]:

$\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(T-T_{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(T-T_{a}\right)+\sigma\left(T^{4}-T^{4}\right)=p} & {r \leq \frac{d_{0}}{2}} \\ {-k \frac{\partial T}{\partial z}+h\left(T-T_{a}\right)+\sigma\left(T^{4}-T^{4}\right)=0} & {r>\frac{d_{0}}{2}}\end{array}\right.$    (6)

And at z=zmax,

$-k \frac{\partial T}{\partial z}=h\left(T-T_{a}\right)+\sigma\left(T^{4}-T_{a}^{4}\right)$    (7)

Eq. (3) alongside with the boundary conditions (Eqns. 4-7) constructs the mathematical model of the heat sink. This set of equations is non-linear 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.

3. Methods and Materials

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. Thermo-physical properties of the heat sink materials used in the study

 Material Thermal conductivity w/m/K Density Kg/m3 Heat capacity J/kg/K Thermal diffusivity m2/s Aluminum type 6061-T6 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 . 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 Lc with passive cooling, facilitate the simplification of the mathematical model by means of the Biot number ($=\frac{h L_{c}}{k}$) . The temperature distribution of the heat sink could be considered uniform for Biot numbers less than 0.1 . For the vertical direction Biotz =0.00015 and for the radial direction Biotr=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 m3 0.08 Thickness of the heat sink - m 0.005 Surrounding (ambient) temperature - K 300 Convective heat transfer coefficient w/m2/K 5 Biot number =hzmax/k vertical Biot number =hd/2/k radial 0.00015 0.0012 Stefan – Boltzmann constant - w/m4/K p=4*Arc tan(1) 3.141592654
4. Lumped Parameter Models

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 convection-radiation 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}=p-h A\left(T-T_{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 Runge-Kutta 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.

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} \theta-c_{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}}}{1-c_{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 hr. The linear form of Eq. (8) is given by:

$\rho c v \frac{d T}{d t}=p-A\left(h+h_{r}\right)\left(T-T_{a}\right)=p-A h\left(1+\xi_{c}\right)\left(T-T_{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(1-e^{-\frac{A h\left(1+\xi_{c}\right)}{\rho c v} t}\right)$                            (15)

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}=p-A\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}{\beta-T}\right)\left(\frac{\beta-T_{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. Parametric Study and Numerical Examples

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 dependent-concentrated solar cell, thermoelectric generator, or high power led. Values of the power input in the range 1-10 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 1-10 w/m2/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 1-5 mm are used for demonstration.

Geometric shape of the cross section of the heat sink: three options for the shape are considered-circular, 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{4-3 \delta}}\right) z_{\mathrm{max}} ; \delta=\left(\frac{a-b}{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/ramanujan-circumference-ellipse/). 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  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).

## 2.png 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).

## 3.png 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).

## 4.png 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

## 5.png Figure 5. Temperature of the heat sink of the convective model and the radiative model vs. time

## 6.png 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).

## 7.png 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).

## 8.png Figure 8. Steady state temperature as a function of ambient temperature

## 9.png 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).

6. Summary and Conclusions

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, 6-9).

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.

Nomenclature
 A area, m2 C specific heat, J. kg-1. K-1 E R T Ta c1, c2 d d0 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, m2. 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 z-1)\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 r-1, 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 z-1)}{\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 r-1, i z)}{\Delta r^{2}}+\frac{\left(T^{n}(i r+1, i z)-T^{n}(i r-1, 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 z-1)}{\Delta z^{2}}$                                (I.6)

Convergence and stability are granted for proper choice of time step .

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 6061-T6

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)

'result-check 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.

References

 Kraus, A.D., Bar-Cohen, A. (1995). Design and analysis of heat sinks, John Wiley and Sons, Inc., New York, (2nd printing 1997).

 Bar-Cohen, A., Kraus, A.D. (1999). Wiley series in thermal management of microelectronics and electronic systems: Incorpera, F. P., Liquid cooling of electronic devices and single-phase convection. John Wiley and Sons, Inc., New York.

 Bar-Cohen, A., Iyengar, M. (2002). Design and optimization of air cooled heat sinks for sustainable development. IEEE CPT Transactions, 25(4): 584-591. https://doi.org/10.1109/TCAPT.2003.809112

 Iyengar, M., Bar-Cohen, A. (2003). Least energy optimization of air cooled heat sinks for sustainable development. IEEE CPT Transactions, 25(4): 16-25. https://doi.org/10.1109/tcapt.2003.811463

 Wang, C.C., Hung, C., Chen, W.H. (2012). Design of heat sink for improving the performance of thermoelectric generator using two-stage optimization. Energy, 39: 236-245. https://doi.org/10.1016/j.energy.2012.01.025

 Lee, S. (1995). Optimum design and selection of heat sinks. Eleventh IEEE SEMI-THERMTM Symposium, 48-54. https://doi.org/10.1109/95.477468

 Aravindh, B.S., Nair, T.R.G. (2008). Heat sink performance analysis through numerical technique. IEEE Symposium (NSSP08), Bangalore, pp. 1-7.

 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): 76-81. https://doi.org/10.1115/1.2429713

 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): 49-54. https://doi.org/10.5755/j01.eee.20.1.6167

 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): 62-72. https://doi.org/10.3390/app7020062

 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 high-power LED bulbs. Journal of Physics: Conference Series, 605: 0120051-0120058. https://doi.org/10.1088/1742-6596/605/1/012005

 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): 61-64. https://doi.org/10.1088/1674.4926/30/4/044011

 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: 164-178. https://doi.org/10.1016/j.solmat.2016.04.016

 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 in-line configurations. The International Journal of Microcircuits and Electronic Packaging, 24(1): 68-76. https://doi.org/10.1088/1674-4926/30/4/044011

 Sun, J., Israeli, T., Reddy, T., Scoles, K., Gordon, J.M., Feuermann, D. (2004). Modeling and experimental evaluation of passive heat sinks for miniature high-flux photovoltaic concentrators. ASME International Solar Energy Conference, 431-440. https://doi.org/10.1115/ISEC2004-65184.

 Smith, G.D. (1985). Numerical Solutions of Partial Differential Equations: Finite Difference Method. 3rd edition, Oxford University Press, NY, 75.

 Smith, G.D. (1985). Numerical Solutions of Partial Differential Equations: Finite Difference Method. 3rd edition, Oxford University Press, NY, 76.

 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.

 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.

 Smith, G.D. (1985). Numerical Solutions of Partial Differential Equations: Finite Difference Method。 3rd edition, Oxford University Press, NY, 43.