Applied Energy 156 (2015) 703–714

Contents lists available at ScienceDirect

Applied Energy
journal homepage: [Link]/locate/apenergy

Multi-objective geometric optimization of a PCM based matrix type

composite heat sink
R. Srikanth, Pavan Nemani, C. Balaji ⇑
Heat Transfer and Thermal Power Laboratory, Department of Mechanical Engineering, Indian Institute of Technology Madras, Chennai 600036, India

h i g h l i g h t s g r a p h i c a l a b s t r a c t

 Experiments conducted on a matrix Numerical model

Matrix type heat sink
type heat sink with n-eicosane as
PCM.
 Full three dimensional simulations
done with Fluent to complement
experiments.
 Simulation driven multi objective
geometric optimization done using
NSGA II.
 Pareto optimal front obtained for the
Artificial neural
problem under consideration.
network
 Optimization increases charging time Pareto optimal solutions
by 5% and reduces discharging time
by 12.5%.

a r t i c l e i n f o a b s t r a c t

Article history: This paper reports the results of a numerical optimization and experimental investigation of phase
Received 13 February 2015 change material (PCM) based composite pin fin matrix heat sink. The main objective of this study is to
Received in revised form 25 June 2015 determine the optimized configuration of the matrix type heat sink that will stretch the operation time
Accepted 15 July 2015
during the heating cycle and minimize the time during the discharging cycle. The PCM used is n-eicosane.
The heat sink is made of aluminium. A constant heat flux of 1.9 kW/m2 is applied at the bottom of the
heat sink. The numerical results were matched up with the experimental results to determine the overall
Keywords:
heat transfer coefficient with the help of commercially available ANSYS FLUENT 14.0 software. For con-
Phase change material
Heat sinks
stant power level and constant volume of the PCM, 40 different geometrical configurations of heat sinks
Electronic cooling were considered and temperature time histories were obtained for both the charging and discharging
Multi objective optimization cycles by using full three dimensional simulations of flow and conjugate heat transfer including phase
change using FLUENT 14.0. The output of these simulations was given as an input to a neural network
and a multi objective optimization was carried out to determine the optimum configuration of the heat
sink, that maximizes the charging period and minimizes the discharging period simultaneously.
Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction

The requirement of compactness has become very critical in the

⇑ Corresponding author.
design of any electronic equipment. Conventional cooling
E-mail address: balaji@[Link] (C. Balaji).

[Link]
0306-2619/Ó 2015 Elsevier Ltd. All rights reserved.
704 R. Srikanth et al. / Applied Energy 156 (2015) 703–714

Nomenclature
L latent heat of PCM, kJ/kg
Greek symbols Tc cooling time, s
a thermal diffusivity, m2/s T liquidus liquidus temperature, K
b thermal expansion coefficient of PCM, K1 Tm melting point, K
c liquid fraction T solidus solidus temperature, K
m kinematic viscosity, m2/s T temperature, K
qs density of Aluminum, kg/m3 Ts time to reach set point temperature, K
t thickness of the heat sink, m
Notations U overall heat transfer coefficient, W/m2 K
Amushy mushy zone constant, m/s v fluid velocity, m/s
cp specific heat at constant pressure, kJ/kg K
Fo Fourier number Abbreviations
Fo h optimum Fourier number – heating Al Aluminum
Fo c optimum Fourier number – cooling ANN Artificial Neural Network
Fon h numerical Fourier number – heating amb ambient
Fon c numerical Fourier number – cooling exp experimental
g acceleration due to gravity, m/s2 GA genetic algorithm
h enthalpy, kJ/kg MRE mean relative error
hs sensible enthalpy, kJ/kg PCM phase change material
hl total latent heat content in body, kg/kg K R2 sum of residuals
H height of the heat sink, m TCE thermal conductivity enhancers
k thermal conductivity, W/m K
ks thermal conductivity of Aluminum, W/m K

techniques are under severe strain to meet the demanding cooling 2. Literature survey
requirements that stem from both the compactness and the
increasing heat fluxes of such electronic devices. The key goal in Pillai and brinkworh [1] considered the use of phase change
the thermal design of any electronic equipment is to keep the material in energy storage systems. They also enumerated the
device temperature within safe operational limits to ensure relia- properties of the phase change materials and summarized on the
bility. Phase change material (PCM) based cooling technique is a selection criteria for a PCM. The authors also emphasized in the
passive type cooling that offers several advantages over conven- design optimization study of a system that uses a PCM. The authors
tional cooling techniques particularly in intermittent duty cycles concluded that PCM offers a convenient means of storing low grade
due to the large latent heat of the PCM. The desirable properties thermal energy within compact systems.
of such a PCM are Voller and Prakash [2] developed an enthalpy-porosity tech-
nique for the numerical modelling for solidification and melting
1. High latent heat of fusion. including the natural convection. The mushy zone was modelled
2. High specific heat capacity. as a porous media and the porosity was set equal to the liquid frac-
3. High thermal conductivity. tion of the corresponding computational grid. Yu and Joshi [3] car-
4. Convenient melting temperature (less than the set point ried out both experimental and numerical studies on natural
temperature). convection, conduction, radiation heat transfer of discrete heat
5. Small volume change during phase change. source using a pin fin heat sink and concluded that enhancement
6. Chemically stability. in heat transfer due to usage of pin fin heat sink was significant.
Akhilesh et al. [4] proposed a thermal design procedure for com-
The major drawback of a PCM based heat sink is however its posite heat sinks (CHS) for maximising the thermal storage and also
low thermal conductivity. As measured the typical conductivity to stretch the operation time of the heat sink under safe operational
of n-eicosane is 0.23 W/m K in the solid state and is much lower limits. Benjamin et al. [5] conducted experimental and numerical
in the liquid state. In a PCM based heat sink the PCM layer located studies on melting of a PCM (n-eicosane) inside a cylinder. The
far away from the heater base can utilize the latent heat only after measured temperatures and the melt front locations of the numer-
the adjoining PCM has fully melted. The low thermal conductivity ical code results were validated against the experimental results.
of PCM, hence delays the full, uniform melting of the PCM, which in They concluded that at a Stefan number of 0.0836 numerical results
turn results in a high temperature gradient within the PCM heat had a good agreement with the experimental data. Nayak et al. [6]
sink. One way of overcoming the high resistance offered by the developed a numerical model for investigating the effectiveness of
PCM to the flow of heat is to add thermal conductivity enhancers using thermal conductivity enhancers to improve the overall ther-
(TCE). The addition of such TCEs ensures lower temperature gradi- mal conductance of the PCM. A transient finite volume method
ent and effective utilization of latent heat without the generation was used to discretise the governing equations.
of hot spots. The effective thermal conductivity of the PCM based Tan et al. [7] experimentally studied the cooling phenomenon
composite heat sinks is also enhanced. The major design goal in in PCM based Personal digital assistants (PDAs) using n-eicosane
a PCM based heat sink would be to stretch the latent heat time dur- as the PCM. The effect of different power levels on both charging
ing charging and to reduce the discharging cycle time. Latent heat and discharging were studied along with the effect of change in
absorption phenomenon with the melting of a particular PCM can orientation. They finally reported that different orientations can
be used to delay the temperature rise of a device subjected to high affect the temperature time histories within the PCM. Martin
heat flux. et al. [8] studied the suitability of PCM for cold storage. This paper
 R. Srikanth et al. / Applied Energy 156 (2015) 703–714 705

discusses the concept of a direct contact PCM water storage for phase change material. The simulations are performed using a val-
peak shaving. Experimental studies showed some drawback of idated mathematical model and an optimal design is arrived at the
using phase change material importantly, its low thermal conduc- end of the simulation.
tivity and volumetric expansion. Yin et al. [9] built an electronic Antunes et al. [21] modeled the sizing and locating the capaci-
cooling experimental system based on thermal adaptation com- tors in electrical radial distribution networks. The aim of the paper
posite material. Here the thermal adaptation composite material was to present an evolutionary algorithm that was particularly sui-
refers to paraffin embedded in graphite. Experimental results ted for thorough non dominated frontier. Among the two objective
showed that the overall heat transfer coefficient is 1.25 times functions that were considered one function is of a technical nature
higher than the traditional cooling system. Furthermore, they con- and the other is of economical nature. The authors reported a set of
cluded that the dosage of composite material has a positive impact non dominated solutions while simultaneously trying to minimize
on the cooling of electronic equipment. installation cost and minimizing the system losses.
Wang and Yang [10] performed three dimensional transient Sanaye and Hajabdollahi [22] performed a multi objective opti-
numerical simulations to investigate the performance of a hybrid mization using non-sorted genetic algorithm (NSGA) to maximize
PCM based heat sink. The number of fins and the heat input values the fin effectiveness and minimize the total annual cost incurred
were varied for each simulation. Available experimental data was to manufacture the fin simultaneously. Six different design vari-
used to validate the numerical model. ables contributing to the conflict of the two objectives were iden-
Hosseinizadeh et al. [11] conducted both experimental and tified, and a correlation between the two objective function and six
numerical studies on the thermal performance of PCM based heat design variables was also developed.
sinks. Studies on the effect of various parameters like number of fins, Ahmed et al. [23] performed a multi objective optimization and
fin height and thickness were studied and an optimum configuration decision making in the selection of a cricket team. Trying to meet
of the fins was proposed. Weng et al. [12] performed an experimen- two objectives, to maximize the batting average and minimize
tal investigation on the performance of heat pipes in conjunction the bowling average, the finite budget was applied as a constraint
with PCM for electronic equipment cooling. The adiabatic section and a set of non dominated solutions was obtained. Sun et al. [24]
of the heat pipe was covered by a storage container with PCM. The proposed a technology that combines phase change materials with
PCM chosen for the current study was tricosane. The authors exper- a cold source for telecommunication base stations. A full prototype
imentally proved that the heat pipe module with PCM could reduce was built and tested in the laboratory. They concluded that the pro-
the heater temperature by 12:3  C as compared to the no PCM case. posed technology helps to reduce a significant amount of energy.
Sertkaya et al. [13] experimentally investigated the effect of ori- Asadi et al. [25] performed a multi objective optimization for build-
entation angle on the natural convection heat transfer of pin finned ing retrofit. Optimization was performed using coupled ANN-GA
heat sink. They concluded that there is a considerable effect on the approach and a set of non dominated solutions was obtained.
heat transfer. Alawadhi and Amon [18] investigated thermal energy manage-
Baby and Balaji [14] conducted experimental investigations on ment issues associated with portable electronic equipment. The
PCM based finned heat sinks using n-eicoanse. The effect of performance of the PCM based thermal control unit was analyzed
different types of fins such as plate fin, pin fin on the melting of for both constant and variable power operations.
PCM was studied and the authors concluded that use of fins helped From the review of literature presented above it is very clear that
in stretching the operational time of the heat sink. the use of phase change material (PCM) based composite heat sinks
Mahmoud et al. [15] experimentally investigated the effect of for electronic cooling is a well established field of research. One can
insert configuration and PCM type on the thermal performance of find that several configurations of heat sinks with PCM has been
the heat sink. The investigation concluded that use of honeycomb investigated experimentally. However, a detailed study of heat
type inserts showed better performance compared to machined transfer in PCM based heat sink due to the presence of pin fins in
metallic inserts. Six different PCMs were used and from the results all the three directions are hard to find in the literature.
it was concluded that the PCM with the lowest melting point was Furthermore, one can find that literature is quite wealthy in opti-
more effective in maintaining the device temperature under the mization studies to maximize the charging period, but one can
safe operating limits. Tari and Mehrtash [16] performed a numer- see little attention has been paid to the discharging cycle of a heat
ical investigation of natural convection heat transfer from inclined sink. This has been considered as a key objective in the current
plate fin heat sinks. The numerical model was validated using ver- manuscript, and this results in the consideration of the
tical plate fins on a vertical base heat sink case. multi-objective optimization. The latter has been accomplished in
Fan et al. [17] experimentally investigated the effect of internal this study using the non sorted genetic algorithm (NSGA II), a state
fins on the performance of the PCM based heat sinks. For the pur- of the art algorithm. During the operation of the heat sink, the
pose of baseline comparison, an unfinned heat sink was taken and charging cycle takes place typically for about 1/5th of the total
it was concluded that use of a PCM could keep the device temper- operating time, while the discharging cycle is 4/5th. This study aims
ature under safe operating limits. Jang et al. [18] performed an to bridge the above mentioned gaps in the literature. Furthermore,
multi objective optimization on a heat sink for LED cooling appli- the objective of this paper is to obtain a set of non dominated solu-
cations. The authors considered both the convection and radiation tions which can maximize the charging period and minimize the
effects. They also conducted experiments to validate the numerical discharging cycle time simultaneously.
model. The two objective functions considered were to minimize The heat sink considered in this study, has a lot of potential in
the thermal resistance and the mass of the heat sink. Finally a the cooling of computer processors. Intel for example uses [26]
set of non dominated solutions were obtained satisfying both the heat sinks with PCMs as a part of its thermal management unit.
objective functions. For the heat sink to be a part of the thermal management sys-
Zhou et al. [19] reviewed the storage of thermal energy using tem the following are some of the key considerations.
phase change material for building applications. The main high-
lights of the paper include the numerical evaluation of thermal  Rotating or moving parts are not preferred due to acoustic
energy storage using PCM in buildings, the criteria for selection constraints.
of PCM, property measurement methods are also discussed.  The heat sink must enable the process to respect the design tem-
Álvaro et al. [20] proposed an innovative design to allow perature limits. For example, for a heat input of 10 W the device
exchange of heat between water and RT 60, which is used as a temperature should not exceed 100 °C after 40 min of operation.
706 R. Srikanth et al. / Applied Energy 156 (2015) 703–714

 The heat sink size needs to be compact in size to accommodate

other cooling equipments like fan, vents and ducts.
 The mass of the heat sink assembly has to be kept minimum so
as to avoid mechanical stress on the PCBs.
 The heat sink assembly should be cost effective.

3. Experimental setup

A pin fin matrix type heat sink (comprising both vertical and hor-
izontal fins) made of aluminium was used in the present study. The
overall dimensions of the heat sink were 80 mm  62 mm  25 mm.
The heat sink cavity has a depth of 20 mm inside which the fins are
present. The fin thickness is 2 mm in all the directions. A wall thick-
ness of 7 mm is maintained in the heat sink. A 2 mm slot is given at
the bottom of heat sink to accommodate the plate heaters. All four
sides of the heat sink are insulated by using a rubber cork and the
top portion of the heat sink is covered with an acrylic. The portion
between the acrylic and heat sink top is sealed by using a silicon rub-
ber gasket. The dimensions of the heater plate are 60 mm  42 mm.
The heater is made of mica sheet and nichrome wire of finite resis-
tance is wound over it.
The pin fin matrix consists of 24 vertical fins, 6 full length hor-
izontal fins in one perpendicular direction and 16 full length hori- Fig. 2. Picture showing the assembled heat sink setup.
zontal fins in other direction, as shown in Fig. 1. The fins acting as
thermal conductivity enhancers (TCE) have a volume fraction of
9%. The volume fraction of 9% TCE was chosen in order to have a
standard comparison with the case of 72 pin fin heat sink which
was earlier performed by Baby and Balaji [14]. The pin fin matrix
was fabricated by the process of Electro Discharge Machining.
The side walls were welded to the matrix by gas welding. A leak
test was performed initially on the heat sink, and when no leak
was observed the experiments were started. The final assembled
heat sink is shown in Fig. 2. A total of 13 calibrated k-type thermo-
couples are used to measure various temperatures in the heat sink.
Calibration was carried out using a constant temperature bath in
the temperature range of 30–70  C. Two thermocouples (T 1 and
T 2 ) are placed at the heat sink base. Four thermocouples (T 3 T 4
T 5 T 6 ) each are placed on each wall. Two thermocouples, each were
placed at 5 mm (T 7 ; T 8 ), 10 mm (T 9 T 10 ), 15 mm (T 11 T 12 ) from the
base inside the PCM. One thermocouple (T 13 ) is placed in still air to
measure the ambient temperature. Fig. 3 shows the positions of

Fig. 3. Picture showing the position of thermocouples in the heat sink.

thermocouples. The thermocouples were fixed to the heat sink

using Araldite epoxy. All the thirteen thermocouples were
connected to AGILENT 34970A data logger. The heating cycle was
carried out until a set point temperature of 45  C is reached.

3.1. Uncertainty analysis

The uncertainty associated with the measurement of funda-

mental quantities like the voltage and current which are the least
counts of the corresponding measurement devices. In so far as the
temperature is concerned, based on the calibration it can be
Fig. 1. Picture showing the heat sink before assembly. inferred that temperature is accurate within ±0.2 °C. However,
 R. Srikanth et al. / Applied Energy 156 (2015) 703–714 707

Table 1 5. Baseline comparison

Results of the uncertainty analysis.

SI no. Quantity measured Uncertainty Unit Baseline comparisons are done in the present study, for the heat
1 Temperature ±0.2 
C sink with PCM against a heat sink without it. For the purpose of
2 Voltage ±0.1 V comparison a new term ‘a’ is defined (enhancement ratio) as the
3 Current ±0.01 A ratio of the time taken by the heat sink to reach a set point temper-
4 Power ±3.2 %
ature with PCM to that without PCM. The baseline results for a
power level of 5 W are shown in Fig. 5. The enhancement ratio
‘a’ for a power level of 5 W is 9.80. As the power level is increased
the derived quantity power requires an uncertainty calculation the enhancement ratio decreases. The enhancement ratios for
which was done and the results are shown in Table 1. other power levels are listed in Table 2.
Uncertainty in power can be calculated as
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
 2  2 6. Material properties
@P @P
rp ¼ rV þ rI
@V @I The PCM material used in the present study is n-eicosane, and
For a power level of 2 W the heat sink material used is aluminium which are same as that
The uncertainty in the power turns out to be used by previous authors [14]. All properties of PCM were taken
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi as constant regardless of the phase except for the specific heat
whose values were measured at 50 temperature points in the
rP ¼  ð0:5  0:1Þ2 þ ðð4  0:01Þ2
range of 30–50  C. In view of this a piecewise linear function for
rp ¼ 0:064 the variation of specific heat with temperature was assumed for
rP
P
¼ 0:032 the purpose of numerical study. The flow is assumed to be laminar
rP ¼ 3:2% and incompressible. The PCM under consideration, n-eicosane is
stable and non-reactive under normal conditions of usage, storage
and transport. The only precaution to be taken while using the PCM
4. Charging and discharging cycles is to avoid contact with strong oxidizing agents, which may cause
instability. The cost of 1 kg of PCM is 230.5 USD.
As discussed earlier, the operation of heat sink takes place in
two cycles, namely charging and discharging cycles. Charging here 7. Numerical model
refers to the process of energizing the heat sink base heater plate
with a constant DC power supply. In the present study, charging The advantage of using simulation to conduct an extensive
is done until the heat sink base reaches a set point temperature parametric study is critical to determine the optimal configuration
of 318.2 K. Discharging here refers to switching off the DC power of the heat sink for maximum performance as the detailed para-
input to the heat sink base. During this process, the heat sink loses metric study through experimentation is time consuming and
heat to the ambient air through natural convection. In the present expensive.
study, experiments were conducted for different power levels
ranging from 3 to 10 W as shown in Fig. 4. As the power input
increases the time taken by the heat sink to reach a particular
set point temperature decreases. As seen from the Fig. 4 for a
power input of 4 W the time taken to reach the set point temper-
ature of 318.2 K is 5620 s, while for a power level of 6 W the time
taken is 2880 s. For a power input of 10 W, the time taken is con-
siderably reduced to 1450 s which is one-fourth of the time taken
for the 4 W case.

Fig. 5. Comparison of temperature–time histories for 5 W power level with PCM to

that of without PCM.

Table 2
Enhancement ratio ‘a’ for 3 power levels.

SI no. Power W a
1 5 9.8
2 8 8.5
Fig. 4. Variation of base temperature with time for 7 different power levels for PCM
3 10 6.9
based matrix type heat sink.
708 R. Srikanth et al. / Applied Energy 156 (2015) 703–714

The modeling of phase-change processes presents a significant

challenge due to complexity and conjugacy of the involved phe-
nomenon [17]. Some important factors that need to be considered
are:

 volume expansion during phase change.

 Convection in liquid phase.
 Motion of solid in the melt due to density differences.

In the current study, the functioning of the heat sink with phase
change material is simulated through the ANSYS FLUENT 14.0 [2].

7.1. Geometry and mesh

Fig. 6 shows the geometry used in obtaining the numerical solu-

tion. The dimensions are the same as those used in experiments.
The model is meshed using tetrahedral meshing on Aluminum
material and hexahedral mesh on [Link] mesh element sizes
are 1 mm and 2 mm on the PCM side and aluminium side respec-
tively. Maximum face size and the maximum size of the grid are
1.2 mm and 2.4 mm respectively. Orthogonal quality of the mesh-
ing is 0.87 and is shown in Fig. 7.

Fig. 7. Figure showing the top view of the meshed model.

7.2. Governing equations

For the Aluminum section, only the conduction equation is Melting is taken care of by incorporating the enthalpy porosity
considered formulation in ANSYS FLUENT 14.0. The enthalpy h can be written
 
@ @ @T as a sum of sensible enthalpy hs and the latent heat hl .
ðqs hÞ ¼ ks þS ð1Þ
@t @xi @xi
h ¼ hs þ hl ð5Þ
In Eq. (1), h is the enthalpy. The section containing PCM is mod-
eled in three dimensional by considering the flow (in molten state) @
to be laminar and compressible flow. ðqhÞ þ $  ðqv hÞ ¼ $  ðk$TÞ ð6Þ
@t
The continuity equation:
where
@q
þ qð$  v Þ ¼ 0 ð2Þ Z T
@t
hs ¼ href þ cp dT ð7Þ
Momentum equations: T ref

@v 1 The liquid fraction, c, can be defined as (source: [2])

þ v  $v ¼  $p þ m$2 v þ g þ S ð3Þ
@t q
c ¼ 0 if T < T solidus ð8Þ
Energy equation:
c ¼ 1 if T > T liquidus ð9Þ
@ T  T solidus
ðqhÞ þ $  ðqv hÞ ¼ $  ðk$TÞ ð4Þ c¼ if T solidus < T < T liquidus ð10Þ
@t T liquidus  T solidus
The latent heat content can now be written in terms of the
latent heat of the material

hl ¼ cL ð11Þ
The enthalpy-porosity technique treats the mushy region (par-
tially solidified region) as a porous medium ([2]). The porosity in
each cell is set equal to the liquid fraction in that cell. In the fully
solidified regions, the porosity is equal to zero, which extinguishes
the velocities in these regions. [2] suggested that the momentum
sink due to the reduced porosity in the mushy zone takes the fol-
lowing form

ð1  cÞ2
S¼ Amush v ð12Þ
ðc3 þ cÞ
Here ‘c’ is a small number (0.001), to prevent the division of the
numerator by zero. When heat is added to a fluid and the fluid den-
sity varies with temperature, a flow can be induced due to the force
of gravity acting on the density variations. Such buoyancy-driven
flows are termed natural-convection (or mixed-convection) flows
Fig. 6. Figure showing the geometry used in modelling. and can be modeled by ANSYS FLUENT 14.0.
 R. Srikanth et al. / Applied Energy 156 (2015) 703–714 709

The modeling is done by involving the Boussinesq approxima- shows the variation of sum of the square of the residuals R2 with
tion. This model treats density as a constant value in all solved heat transfer coefficient for the heating cycle and a similar plot
equations, except for the buoyancy term in the momentum was obtained for the cooling cycle also. However to develop the
equation numerical model, five thermocouple measurements which include
the temperature reading inside the PCM was also considered to
ðq  q0 Þg q0 bðT  T 0 Þg ð13Þ
estimate the overall heat transfer coefficient (U). The value of U
The Boussinesq approximation is given as follows is determined by minimizing the sum of the square of the error
in the temperatures in least square sense. Thus the final value of
q ¼ q0 ð1  bTÞ ð14Þ
U has the convection effects included in it, which is the base for
This eliminates the q from the buoyancy approximation developing the numerical model. Therefore the heat transfer coef-
The change in material properties of the material with melting ficients are chosen to be 1.5 and 8 W/m2 K for the heating and cool-
are modeled using a piecewise linear model for the specific heat ing cycles respectively. The heat transfer coefficient for the cooling
and thermal conductivity. SIMPLE algorithm is used with pressure cycle is higher because of the average temperature during the cool-
velocity coupling. PRESTO scheme is adopted for discretizing of ing cycle is higher than that of the heating cycle leading to higher
pressure and SIMPLE algorithm with the second order upwind natural convection coefficients.
scheme is used. The convergence is set to 1  106 ; 1  103 and Initial conditions At t = 0 (was set in accordance to the
1  106 for the continuity, momentum and energy equations experiment)
respectively. An energy balance accuracy of 4:332  109 was
T ¼ T Al ¼ T F ¼ T S ¼ 305:5 K ð15Þ
achieved after convergence.
For the top surface of the heat sink
7.3. Grid independence studies
@T @T Al @T F @T S
¼ ¼ ¼ ¼0 ð16Þ
The heat sink is meshed with tetrahedral meshes on the @y @y @y @y
Aluminum material and with hexahedral mesh on the PCM. To
determine the number of elements that needs to be taken into for all the other outer surfaces
account for simulation, a grid independence study was conducted.
Fig. 8 shows the results of the grid independence studies and it is @T Al
¼0 ð17Þ
seen that a mesh with about 284,000 nodes is adequate and this @y
was used for all subsequent computations.
For the surface in contact with the heater

7.4. Boundary conditions @T Al

kAl ¼q ð18Þ
@y
The heat sink is subjected to two cycles viz. heating cycle and
the cooling cycle. During the heating cycle, the heat sink is sub-
jected to a constant power input of 5 W.
Heat losses occur around the heat sink owing to natural convec- 7.5. Validation
tion and these need to be quantified. Although the heat transfer
coefficient varies with the surface temperature, an overall average The numerical results for the temperature time history of the
heat transfer coefficient is determined by comparing the tempera- base are compared with the experimental result and the results
ture time histories simulated for various values of U with respect are shown in Fig. 10. The close match between the two is evident.
to the experimental temperature time history. For this, the heating As already mentioned, the heat sink is heated till it reaches a set
cycle and the cooling cycle are treated separately and the overall point temperature of 45  C. By then all the PCM would have molten
heat transfer coefficient is varied in steps of 0.5 W/m2 K. Fig. 9 completely as the melting temperature of the PCM is 36:5  C.

Fig. 9. Variation of sum of residuals, S with heat transfer coefficient, h for heating
Fig. 8. Variation of base temperature with time for different meshes. cycle.
710 R. Srikanth et al. / Applied Energy 156 (2015) 703–714

Table 3
Geometrical models generated using LHS method.

Model No. Length mm Width mm Height mm

1 62.9 32.0 31.5
2 65.3 47.2 20.6
3 60.5 37.6 27.9
4 59.3 33.6 31.8
5 73.1 50.4 17.2
6 61.1 39.2 26.5
7 67.1 39.2 30.3
8 64.1 44.0 22.5
9 55.1 44.8 22.1
10 58.7 48.8 22.1
11 56.9 48.0 23.2
12 53.9 43.2 27.2
13 66.5 52.0 18.3
14 56.3 49.6 22.7
15 64.7 52.8 18.5
16 71.3 42.4 21.0
17 68.3 34.4 27.0
18 52.1 45.6 26.7
19 70.1 40.8 22.2
Fig. 10. Comparison of numerical and experimental base temperature for PCM 20 71.9 53.6 16.4
based matrix type composite heat sink.

8. Optimization

The goal of this investigation is to determine the optimal

configuration of the heat sink with two objectives viz maximizing
the heating time and minimizing the cooling time, given the con-
straint that the volume of the phase change material is constant.
The following methodologies are adopted to carry out the
optimization.

1. Generate geometries with random dimensions.

2. Perform numerical simulations for each geometry.
3. Train a neural network to find a input–output correlation.
4. Input a population of data and using multi optimization tool
and obtain the optimal outcomes in the form of a Pareto plot.

8.1. Latin Hypercube Sampling

Fig. 11. Variation of base temperature with time for 4 models.
The geometries are generated using the Latin Hypercube
Sampling (LHS) technique. LHS is a statistical method for generat-
ing a sample of the plausible collection of parameter values from a
multidimensional distribution. A sampling function of N variables charging and discharging cycle. Hence, there appears to be a excel-
is divided into M intervals. M sample points are then placed to lent scope for multi objective optimization.
satisfy the Latin Hypercube requirements. In the present study,
the volume of the PCM cavity is a constant and so the independent 8.2. Artificial Neural Network
parameters are length and height of the heat sink (N = 2) with the
width getting fixed automatically. M corresponds to the number of An Artificial Neural Network (ANN) is a non linear regression
geometries generated which is 20 in this study. It is instructive to tool inspired by the biological neuron. It consists of an input layer,
note the symmetry of the problem in terms of the length and hidden layer and an output layer. The hidden layer consists of
width. This in effect gives us 40 inputs to train the neural network. nodes which assigns weights to their inputs and add some biases.
As the length and height are varied the spacing between the fins These are decided while training the neural network for the given
adjust accordingly keeping the percentage of TCE constant as 9% input–output data. The architecture of ANN employed in the pre-
The dimensions of the geometries thus generated are shown in sent study is shown in Fig. 12.
Table 3. The inputs to the ANN are the length (l) and height (h) of the
For each of these geometries, simulations are done for the heat- PCM and the outputs are the heating time (T h ) and the cooling time
ing and cooling cycles to study the fluid flow and heat transfer (T c ) represented non-dimensionally as the Fourier number (Fo)
behavior. The variation of base temperature for select 4 models is which is a function of the geometry of the heat sink viz. length
shown in Fig. 11. From the figures it can be observed that the case (l), height (h) power level (Q) and the total volume of the PCM.
18 takes a much longer time to completely solidify as compared to The Fourier number is defined as
the remaining cases. Further, the case 4 reaches the set point tem- t
Fo ¼ aPCM ð19Þ
perature 150 s later compared to the other cases. From the 40 mod- H2
els that were generated and monitoring their corresponding Fo ¼ f ðl; b; Q ; V t Þ ð20Þ
outputs, there seems to be a considerable diversity spread in both
 R. Srikanth et al. / Applied Energy 156 (2015) 703–714 711

Table 4
Results of the Neuron independent study.

Fig. 12. Architecture of the neural network employed in the present study.

The optimization is carried out for the PCM volume (V t ) to be

63,360 mm3 for a constant power supply of 5 W and TCE to be at
9%. Further, the bounds for the dimensions of the PCM are limited X
m

as follows SðbÞ ¼ ½yi  f ðxi ; bÞ ð29Þ

i¼1
30 mm 6 l 6 80 mm ð21Þ
Of the 40 data points, 80% data points are used to train the network
15 mm 6 h 6 30 mm ð22Þ and the remaining data points are used to validate the network.

The input variables to the ANN are non-dimensionalized as l

and h which are defined as 9. Genetic algorithms

 l Genetic algorithms belong to the class of evolutionary algo-

l ¼ ð23Þ
lmax rithms which generate solutions to optimization problems that
 h are inspired by natural phenomenon like natural selection and
h ¼ ð24Þ
hmax inheritance. Genetic algorithms show great ability to determine
the global optimum in a domain consisting of many local optima.
where lmax ¼ 80 mm and hmax ¼ 30 mm. Thus, the This is done by initiating the search at multiple points in the search
non-dimensionalized values are in the following range space and subsequent operations such as mutation, crossover and
0:375 6 l 6 1

ð25Þ selection are implemented till the optimum is obtained, so that a
 robust and extensive search of the parameter space is done.
0:5 6 h 6 1 ð26Þ
A neuron independence study was conducted to determine the 9.1. Multi objective optimization
optimum number of neurons to be used in the hidden layer.
Performance parameters such as the Root mean square error In the present study two objectives, namely maximizing the
(RMS) and mean relative error (MRE) defined below, were calcu- heating time and minimizing the cooling time are considered.
lated for networks constructed with varying nodes in the hidden
Objectiv e function f 1 ¼ maximize charging time
layer. The parameters are defined as
Objectiv e function f 2 ¼ minimize discharging time
PN
 tnum;i j
i¼1 jt ANN;i
MRE ¼ ð27Þ Subject to constraints
N
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
PN 2 Cav ity v olume ¼ 63; 360 mm3
i¼1 ðt ANN;i  t num;i Þ
RMS ¼ ð28Þ TCE v olume ¼ 6336 mm3
N
Table 4 shows the variation of MRE and RMS with number of In contrast to a single objective optimization, it may not be pos-
neurons in the hidden layer. From this study, it was found that a sible to have a unique solution satisfying both the objectives. So a
total of 9 neurons in the hidden layer is adequate. trade off is necessary and this leads us to the concept of Pareto
A two layer feed forward back propagation network is used and optimal solutions. A solution is called Pareto optimal or non dom-
the network is trained using the Levenberg Marquadart algorithm inated if none of the objectives can be improved without degrading
in MATLAB software. The Levenberg Marquadart algorithm is used other objectives. In mathematical terms, a multi objective opti-
for a least squares curve fitting problem. For a given set of m mization problem can be formulated as
empirical data points, input (x) and output (y) variables, the minðf 1 ðxÞ; f 2 ðxÞ; . . . ::; f k ðxÞÞ ð30Þ
optimized value b is derived such that the sum of squares of the
g i ðx 6 0ðj ¼ 1; 2; . . . :; mÞÞ ð31Þ
deviations (SðbÞ) is minimized.
712 R. Srikanth et al. / Applied Energy 156 (2015) 703–714

where x 2 X having k objectives and m constraints (by function g). It is very clear from the figure that it is impossible to move from
A solution x1 is said to Pareto dominate another solution x2 if one point to another without degrading any one of the objective
function. For instance,data points belonging to the top most cluster
f i ðx1 Þ 6 f i ðx2 Þ 8i 2 ð1; 2; . . . kÞ ð32Þ A in the front satisfy the objective function 1 more compared to the
1 2 objective 2. In contrast the data points belonging to cluster D
f j ðx Þ < f j ðx Þ for at least one j 2 ð1; 2; . . . kÞ ð33Þ
satisfy the objective 2 more compared to the objective 1. If more
A solution which is not dominated by any other solution in the weightage is given to objective 1 points higher on the Pareto front
solution space is known as Pareto optimal solution. A set of all need to be chosen and for a higher weightage of objective 2 will be
these non dominated solutions is called the Pareto optimal set vice versa. The multi objective optimization yielded a set of non
and the corresponding function values in the objective space are dominated solutions which could simultaneously increase the
called as the Pareto front. charging time and decrease the discharging cycle time, which
In order to do this, a non sorted genetic algorithm (NSGA) [23] is could be treated as a first of a kind investigation on a heat sink with
used. NSGA is an extension of genetic algorithm for multi objective PCM for electronic cooling. Another important implication of
optimization. The objective of the NSGA is to improve the adaptive Fig. 13 is that 80% of the points generated on the Pareto front
fit of a population of candidate solutions to a Pareto front con- corresponds to the geometries of high surface area compared to
strained by a set of objective functions. This is different from the the initial geometries. Hence, the increase in surface area enhances
GA only in the way the selection operator works. The population the heat transfer during the discharging cycle and cooling time is
is sorted into a hierarchy of sub-populations based on the ordering minimized. Furthermore the geometries corresponding to the
of Pareto dominance. Similarity between members of each points on the Pareto front have optimum spacing between the fins,
sub-group is evaluated on the Pareto front, and the resulting which helps in full utilization of the latent heat of the PCM, thereby
groups and similarity measures are used to promote a diverse front stretching the time of operation of the device.
of non-dominated solutions.
In the current study, objective 1 is the time required for heating
9.2. Validation of optima
and objective 2 is the time required for cooling. The goal is to max-
imize objective 1 and minimize objective 2 simultaneously. In order to have more faith in the optimization methodology
MATLAB 2013b is used to implement the NSGA for multi objective
proposed here, the numerically generated optima using the
optimization. Since this solves for a minimization problem, nega- ANN-GA was validated by performing numerical simulations. The
tive of objective 1 is used. Over 400 data points were generated
optimized configurations were verified using the validated
to determine the Pareto front. Fig. 13 shows the Pareto front numerical model that was developed using ANSYS FLUENT 14.0.
obtained for the multi objective problem. The figure shows a set
The corresponding length, breadth and height were modeled and
of non dominated solutions, each of which satisfies both the objec- the estimated overall heat transfer coefficient was applied. The
tives. 32 distinct points were obtained on the Pareto front. Four
time taken to reach the set point temperature in heating was
random Pareto optimal points were chosen from each of the clus- 3995 s, which was very close to the Fourier number of 2.4
ters A, B, C and D as shown in Fig. 13. The clustering of the data set
predicted by hybrid ANN-GA. The results predicted by the
was performed using the well known ‘‘K-Means clustering’’ algo- optimization show a good agreement as shown in Fig. 14.
rithm. The algorithm breaks the data set into k different clusters.
‘k’ random observations are chosen from the data set and are
9.3. Fluid flow and heat transfer characteristics of optimal
assigned as the seeds initially. All other observations are combined
configuration
with the seeds based on the proximity of those observations from
the seeds. The algorithm aims to minimize the distance between
The temperature contours at the mid plane along the length for
the observation data point and the seed point. After the set of clus-
the optimized case l = 51.46 mm, b = 46.38 mm, h = 26.28 mm (one
ters are formed the center point or the mean point of each cluster is
among the non dominated solutions) for a heat input of 5 W are
assigned as the seeds in the next iteration. The process is continued
shown in Fig. 15. Heat conduction takes place from the heat sink
until the mean in the previous iteration and the current iteration
body to the nearby PCM. At a particular time step high
does not show much variation.

Fig. 14. Figure showing the comparison of t Fourier number (ANN-GA) to the
Fig. 13. Pareto front. Fourier number (numerical) in heating cycle.
 R. Srikanth et al. / Applied Energy 156 (2015) 703–714 713

Fig. 15. Isotherms for optimized heat sink at time (i) 1000 s, (ii) 2000 s, (iii) 3000 s.

temperatures are formed near the aluminium surface. As observed

from the isotherms, the layer of the PCM adjacent to the heater
base gets heated up first and the temperature of that layer shoots
up. Once the layer gets fully melted (reaches a temperature more
than the melting point), the adjacent layer is heated up and the
heat flow progresses. The propagation of heat is much faster in this
matrix type heat sink due to the presence of both horizontal and
vertical fins. From the liquid fraction contours it is evident that
the PCM near the heat sink base and some portions near the wall
has completely melted at the end of 1000 s as shown in Fig. 16.
At the end of 2000 s, the conduction of heat in both horizontal
and vertical direction aids the 50% of the PCM to melt. This is an
excellent sign which indicates the utilization of PCM latent heat
to a larger extent.
At the end of 3000 s, complete melting of PCM is observed. The
disadvantage of having a matrix type heat sink can thus be inferred
from the fact that the entire PCM has melted at the end of 3200 s,
which limits the latent heat time of the entire heat sink to 3200 s.
Fig. 16. Liquid fraction contours for optimized heat sink at time (i) 1000 s, (ii)
The merit of PCM based heat sink is largely based on the latent heat 2000 s, (iii) 3000 s.
period, which for a matrix type heat sink is limited due to the dif-
fusion of heat in all the three directions. However the operation
time in this study is defined as the time taken to reach a set point During melting the effect of convection is present only when
of 318.2 K, which for a matrix type heat with optimized configura- the Stefan number is greater than 1. For a constant heat flux at
tion turns out to be more than 4100 s. the base, the Stefan number can be calculated as [27].
In this study, the fin cross section is held constant, whereas the Stefan number
overall length, breath and height is varied. This indirectly is an
q0 C p H
optimization of the spacing between the fins in X, Y and Z direction. Ste ¼
kp L
In Fig. 6 L, B and H are the overall dimensions. Sb, SL are the spacing
between the fins in the X, Y directions respectively. Therefore the Stefan number in this study is chosen more than
The spacing between the fins plays an important role in the 1, by choosing the appropriate heat input at the base (5 W) to
melting and solidification phenomenon as the mode of heat trans- account for the effect of natural convection on melting process.
fer during melting is both conduction and convection, while in The Stefan number denotes the amount of superheating experi-
solidification the only mode of heat transfer is conduction. If only enced by the liquid phase of the PCM.
melting is considered the spacing should be optimum enough to During the solidification, it is better to have the Aluminum sur-
allow sufficient area for the convection currents, which originate face close to each other so that more conduction takes place.
from the heated surface(base, walls, fin surface) and which go on As soon as the heat flux at the bottom is made zero, the liquid
and erodes the unmelted solid phase nearby and enhance melting PCM near the wall solidifies and forms a layer for solid PCM (con-
of PCM. duction resistance). Due to low thermal conductivity of PCM, the
714 R. Srikanth et al. / Applied Energy 156 (2015) 703–714

heat from the inside has to be conducted through this low thermal References
conductivity PCM layer and then reach the aluminium walls. As the
time increases, resistance to the conduction is increased by the [1] Pillai K, Brinkworth B. The storage of low grade thermal energy using phase
change materials. Appl Energy 1976;2(3):205–16.
subsequent formation of solid PCM layers. In view of the above, [2] Voller V, Prakash C. A fixed grid numerical modelling methodology for
it is clear that, the solidification process needs less spacing convection-diffusion mushy region phase-change problems. Int J Heat Mass
between fins for effective conduction of heat, as the PCM layer Transfer 1987;30(8):1709–19.
[3] Yu E, Joshi Y. Heat transfer enhancement from enclosed discrete components
formed has a very low thermal conductivity. So there is a high using pin–fin heat sinks. Int J Heat Mass Transfer 2002;45(25):4957–66.
scope for multi objective optimization of the spacing between [4] Akhilesh R, Narasimhan A, Balaji C. Method to improve geometry for heat
the fins. Now when both melting and solidification is considered transfer enhancement in PCM composite heat sinks. Int J Heat Mass Transfer
2005;48(13):2759–70.
the GA arrives at a set of non dominated solutions which will be
[5] Jones BJ, Sun D, Krishnan S, Garimella SV. Experimental and numerical study of
the optimum spacing between the fins, that satisfy both the melting in a cylinder. Int J Heat Mass Transfer 2006;49(15-16):2724–38.
objectives. [6] Nayak K, Saha S, Srinivasan K, Dutta P. A numerical model for heat sinks with
phase change materials and thermal conductivity enhancers. Int J Heat Mass
Transfer 2006;49(11-12):1833–44.
10. Conclusions and outlook [7] Fok S, Shen W, Tan F. Cooling of portable hand-held electronic devices using
phase change materials in finned heat sinks. Int J Therm Sci
2010;49(1):109–17.
Experiments were conducted on a phase change material (PCM) [8] Martin V, He B, Setterwall F. Direct contact PCM water cold storage. Appl
based matrix type heat sink. A numerical model for the exactly the Energy 2010;87(8):2652–9.
same geometry was developed using commercially available [9] Yin H, Gao X, Ding J, Zhang Z, Fang Y. Thermal management of electronic
components with thermal adaptation composite material. Appl Energy
ANSYS FLUENT 14.0. Numerical results were matched up with
2010;87(12):3784–91.
experimental results for the case of uniform heat input of 5 W to [10] Wang Y-H, Yang Y-T. Three-dimensional transient cooling simulations of a
estimate the overall heat transfer coefficient which was then used portable electronic device using PCM (phase change materials) in multi-fin
heat sink. Energy 2011;36(8):5214–24. PRES 201.
in the optimization studies.
[11] Hosseinizadeh S, Tan F, Moosania S. Experimental and numerical studies on
Subsequent to this, matched simulations of fluid flow and heat performance of PCM-based heat sink with different configurations of internal
transfer were done for 40 different geometric configurations of the fins. Appl Therm Eng 2011;31(17-18):3827–38. 2010 Special Issu.
heat sinks (30 mm 6 l 6 80 mm; 15 mm 6 h 6 30 mm). [12] Weng Y-C, Cho H-P, Chang C-C, Chen S-L. Heat pipe with {PCM} for electronic
cooling. Appl Energy 2011;88(5):1825–33.
The important conclusions from the present study are as [13] Sertkaya AA, Bilir A, Kargici S. Experimental investigation of the effects of
follows. orientation angle on heat transfer performance of pin-finned surfaces in
natural convection. Energy 2011;36(3):1513–7.
[14] Baby R, Balaji C. Experimental investigations on phase change material based
1. For a PCM based matrix type heat sink, it is possible to use finned heat sinks for electronic equipment cooling. Int J Heat Mass Transfer
NSGA II multi-objective optimization approach to optimize 2012;55(5-6):1642–9.
the heat sink dimensions that can stretch the operation time [15] Mahmoud S, Tang A, Toh C, AL-Dadah R, Soo SL. Experimental investigation of
inserts configurations and PCM type on the thermal performance of PCM based
of the heat sink during the charging cycle by 5% and reduce heat sinks. Appl Energy 2013;112(0):1349–56.
the discharging cycle time by 12.5% simultaneously. [16] Tari I, Mehrtash M. Natural convection heat transfer from inclined plate-fin
2. Using NSGA II multi-objective optimization algorithm, a set of heat sinks. Int J Heat Mass Transfer 2013;56(1-2):574–93.
[17] Fan L-W, Xiao Y-Q, Zeng Y, Fang X, Wang X, Xu X, Yu Z-T, Hong R-H, Hu Y-C,
non-dominated solutions were obtained instead a single opti- Cen K-F. Effects of melting temperature and the presence of internal fins on the
mum solution. The shape of the Pareto curve indicated the performance of a phase change material (PCM)-based heat sink. Int J Therm Sci
trade-off between the two objectives. Upon the analysis of the 2013;70(0):114–26.
[18] Jang D, Yook S-J, Lee K-S. Optimum design of a radial heat sink with a fin-
results, the physical characteristics of the phenomenon are well
height profile for high-power led lighting applications. Appl Energy
portrayed and also helps us to understand the simultaneous 2014;116(0):260–8.
influence of the input variables (l, b, h) on the performance of [19] Zhou D, Zhao C, Tian Y. Review on thermal energy storage with phase change
the heat sink. The cases corresponding to the optimal solution materials (PCMs) in building applications. Appl Energy 2012;92(0):593–605.
[20] Campos-Celador Ãlvaro, Diarce G, Zubiaga JT, Bandos TV, Garcia-Romero AM,
from each cluster was validated using the numerical model Lapez L. Design of a finned plate latent heat thermal energy storage system for
developed and excellent agreement was found. domestic applications. Energy Procedia 2014;48(0):300–8. proceedings of the
3. This methodology adopted results in a set of non dominated 2nd International Conference on Solar Heating and Cooling for Buildings and
Industry (SHC 2013).
solutions. The algorithm can be extended to a wide class of [21] Antunes CH, Pires DF, Barrico C, Gomes Ãlvaro, Martins AG. A multi-objective
problems where optimization is required to satisfy two objec- evolutionary algorithm for reactive power compensation in distribution
tives which are equally important. networks. Appl Energy 2009;86:977–84.
[22] Sanaye S, Hajabdollahi H. Thermal-economic multi-objective optimization of
plate fin heat exchanger using genetic algorithm. Appl Energy
The results from the present investigation on multi objective 2010;87(6):1893–902.
geometric optimization of a matrix type heat sink with PCM can [23] Ahmed F, Deb K, Jindal A. Multi-objective optimization and decision making
approaches to cricket team selection. Appl Soft Comput 2013;13(1):402–14.
[24] Sun X, Zhang Q, Medina MA, Liu Y, Liao S. A study on the use of phase change
1. Serve as a benchmark for future numerical investigations. materials (PCMs) in combination with a natural cold source for space cooling
2. Pave way for more focused research on optimization studies in telecommunications base stations (TBSs) in china. Appl Energy
2014;117(0):95–103.
with simultaneous consideration of charging and discharging
[25] Asadi E, da Silva MG, Antunes CH, Dias L, Glicksman L. Multi-objective
cycles. optimization for building retrofit: a model using genetic algorithm and
artificial neural network and an application. Energy Build 2014;81(0):444–56.
Future studies can consider flow visualization with a viewing [26] Electronic package technology development. Intel Technology Journal,
[Link]. 2005;(0):1–9.
window in the setup, benchmarking simulations and articulating [27] Shatikian V, Ziskind G, Letan R. Numerical investigation of a PCM-based heat
the reasons for the optimum performance from a fundamental sink with internal fins: constant heat flux. Int J Heat Mass Transfer
viewpoint. 2008;51(0):1488–93.