International Journal of Heat and Mass Transfer 55 (2012) 5514–5525

Contents lists available at SciVerse ScienceDirect

International Journal of Heat and Mass Transfer

journal homepage: www.elsevier.com/locate/ijhmt

Optimum tip gap and orientation of multi-piezofan for heat transfer enhancement
of finned heat sink in microelectronic cooling
M.K. Abdullah a,⇑, N.C. Ismail b, M. Abdul Mujeebu c, M.Z. Abdullah b, K.A. Ahmad d, Muhamad Husaini b,
M.N.A. Hamid a
a
Mechanical Section, Universiti Kuala Lumpur, Malaysian Spanish Institute, 09000 Kulim Hi-Tech Park, Kedah, Malaysia
b
Aerodynamic and Advanced Cooling Laboratory, School of Mechanical and Aerospace Engineering, Engineering Campus, Universiti Sains Malaysia, 14300 Nibong Tebal,
Penang, Malaysia
c
Department of Mechanical Engineering, Anjuman Institute of Technology and Management (Visvesvaraya Technological University, Belgaum), 58320 Bhatkal, Karnataka, India
d
Department of Aerospace Engineering, Faculty of Engineering, Universiti Putra Malaysia, 43400 Serdang, Selangor, Malaysia

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

Article history: Piezoelectric fans can be manipulated to generate airflow for cooling microelectronic devices. Their out-
Received 28 September 2011 standing features include noise-free operation, low power consumption and suitability for confined
Received in revised form 6 May 2012 spaces. This paper presents experimental optimization of tip gap and orientation angle of three piezoelec-
Accepted 6 May 2012
tric fans (multi-piezofan) to maximize the heat removal performance of finned heat sink for microelec-
Available online 19 June 2012
tronic cooling. Design of experiments (DOE) approach is used for the optimization, and a three
dimensional simulation using FLUENT 6.3.2 is carried out to better understand the flow induced by the
Keywords:
multi-piezofan and the resulting heat transfer from the heat sink surface. For the optimization, the Cen-
Finned heat sink
Multi-piezofan
tral Composite Design (CCD) of response surface methodology (RSM) is exploited from the Design Expert
Tip gap software. In the numerical model, the flow induced by the piezofan is treated as incompressible and tur-
Optimization bulent; the turbulence is taken care by the shear stress transport (SST) k–x model. The experimental
Temperature drop results are found to be in good agreement with the predictions. Out of 13 experimental trials determined
Heat transfer coefficient by CCD, the optimum tip gap and fan orientation are found to be d = 0.17 and 90° respectively. At this
condition, an enhancement in convective heat transfer coefficient exceeding 88% is achieved, compared
to natural convection.
Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction with a temperature drop of 17 °C. Ihara and Watanabe [4] investi-
gated the flow around the ends of oscillating flexible cantilevers.
Piezoelectric fan (named henceforth as piezofan) is a cantilever The experimental flow profiles were matched with the simulation
beam bonded with a piezoelectric material near its clamped end. flow fields developed by discrete vortex method. Acıkalın et al. [5]
The mechanism of alternate expansion and contraction generates compared analytical, computational and experimental flow pro-
shear forces within the beam, whenever an alternative input signal files for the baffled fan, and obtained a close match among them.
is applied to the piezoelectric material. These forces produce oscil- Advanced flow measurements were carried out by Abdullah et al.
lations at the free end of the cantilever beam which creates motion [6,7] at different piezofan heights by using particle image veloci-
in the surrounding fluid. As piezofans facilitate noise-free opera- metry (PIV) system.
tion with low power consumption, and require less space, they Schmidt [8] focused on the local and average transfer coeffi-
can be adapted to operate at frequencies which are too quiet to cients on a vertical surface by using two piezofans in out-of-phase,
the human ear and can be built and modified to meet various geo- and showed that the transfer coefficients were affected by the tip
metric constraints for many applications. Owing to the excellent gap, and the fan spacing. In a similar study using single piezofan,
thermal management capability, recently piezofans gained signifi- Acikalin et al. [9] obtained significant localized cooling and more
cant attention for electronic cooling applications [1]. than 100% enhancement in convective heat transfer coefficient rel-
Piezofan was first emerged in the early seventies, and was rein- ative to natural convection. In their extension work [10], it was
troduced by Toda [2,3] who found that the piezofans were able to shown that for a given power consumption piezofans were signif-
cool either side of a power transistor panel of a television receiver, icantly better than axial fans, and compared to natural convection
heat sink the piezofans used lesser volume. Yoo et al. [11]
⇑ Corresponding author. attempted the possibility of replacing rotary fan by piezofan.
E-mail address: mkhalil@eng.usm.my (M.K. Abdullah). Different vibrating metal plates were analyzed theoretically and

0017-9310/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved.
http://dx.doi.org/10.1016/j.ijheatmasstransfer.2012.05.024
 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525 5515

List of symbols

A movement area (m2) tp piezofan thickness (mm)

Af frequency dependent amplitude ~
u velocity vector
Dx cross-diffusion term u velocity (m/s)
F external body force (N) ug grid velocity
G gap between piezofan and heat sink (mm) V control volume
e
G generation of turbulence kinetic energy w beam position (mm)
k
ex
G generation of x w_ beam velocity (rad/s)
g gravitational acceleration (m/s2) wp width of piezofan (mm)
Hf fin height (mm) x, y, z space coordinates
hhs heat convection coefficient (W/m2 K) Y dissipation coefficient
kt turbulent kinetic energy (m2/s2)
lp length of piezofan (mm) Greek symbols
lu length of un-patch piezofan (mm) a angle
P static pressure (N/m2) b eigen value
q_ h heater heat generation (W/m3) d normalized tip gap
R2 correlation coefficient e turbulent kinetic energy dissipation rate (m2/s3)
S source term q fluid density (kg/m3)
SD standard deviation s viscous stress tensor (N/m2)
T temperature (K) C diffusion coefficient
Tj junction temperature (K) U general scalar for transport equation
Tn=1,2.. temperature at a n-point (K) x ratio of e to k (s1)
Ts-exp experiment heat sink surface temperature (K) xb beam angular velocity (rad/s)
T1 ambient temperature (K)
t time (s)
th heater thickness (mm)

showed that piezofans were capable of generating jet flow for cool- a constant heat flux of 9166.7 W/m2 is applied. One side of the hea-
ing purpose. Hosaka and Itao [12] and Basak and Raman [13] re- ter is embedded on a wood platform which also serves as thermal
ported that the vibration characteristics of a vibrating cantilever insulator, while the other side is pasted to the heat sink bottom sur-
were altered by the presence of a second oscillating beam depend- face, by using RSÒ heat sink compound which has high thermal con-
ing on the vibration amplitude as well as the pitch, and the phase ductivity. The target being cooled is an aluminum heat sink
difference between the neighboring cantilevers. These findings (Fig. 2(a)), which has base plate of size 80 mm  60 mm  1 mm,
were also confirmed experimentally by Kimber and co-workers and is provided with four fins of 1 mm thickness and 30 mm height
[14,15] who studied the fluidic coupling between piezofans. (total footprint of the heat sink is 80 mm  60 mm  31 mm). Eight
While most of the researchers focused on the use of single K-type thermocouples are used to monitor temperatures at critical
piezofan, only few works [4,8,12–16] were reported on multi- locations as shown in Fig. 2(b) and (c). An additional thermocouple
piezofan which has important practical applications. Few research- is used to monitor the ambient air temperature in the enclosure.
ers were interested in the combination of piezofan and finned heat The temperatures are displayed and recorded by the computer by
sink [17–19], which was proved to be relatively excellent. However means of data logger, for each minute during the experiment.
the combination of multi-piezofan with finned heat sink is yet to The multi-piezofan consists of three piezofans that can be ori-
be explored. Moreover, a three dimensional (3D) numerical analy- ented at angles (a) of 0° to 90°; Table 1 shows the specifications
sis on the heat transfer performance of piezofans is also lacking.
Accordingly, the present study is focused on experimental analysis
and optimization, and 3D modeling of multi-piezofan arranged in
conjunction with finned heat sink. The coupled effect of tip gap
(d) and fan orientation angle (a) on the heat transfer performance
of the heat sink is studied. As an excellent tool for experimental de-
sign and optimization, the design of experiment (DOE) approach is
employed for the estimation of the number of experimental trials,
and subsequent optimization of d and a with the objective of min-
imizing the temperature (maximize the temperature drop) of the
heat sink surface; the procedure was similar to that of Elnaggar
et al. [20]. Based on the best and worst cases of d and a, 3D simu-
lations are performed by using FLUENT 6.3.2, and the results are
compared with the experimental findings.

2. Experimental setup and procedure

Fig. 1 shows the experimental apparatus for the temperature

measurement. The test section of size 800 mm  120 mm 
250 mm is made from clear glass. The electronic package is
mimicked as a heater (size 30 mm  30 mm  3 mm) into which Fig. 1. The experimental setup to study the temperature.
5516 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525

Table 1
Specifications of the piezoelectric fan.

Specification Value
Material Stainless steel
Fan size (mm) 47 (lP)  12 (wP)  0.4 (tP)
Resonant frequency (Hz) 111
Power consumption of fan and circuit (mW) 42
Fan weight (gm) 2.0

Fig. 3. (a) Fan height and (b) Piezofan gap height and orientation angle, and fin
height.

be used, which is defined as the ratio of gap height (G) to the fin

height (Hf); a, Hf and G are illustrated in Fig. 3.
A laser displacement system (KEYENCE LK-G152) was used to
measure the deflection produced by the piezofan; Fig. 4 shows
the experimental setup. The laser beam is pointed on the tip of
the vertically mounted piezofan. When the piezofan swings left
and right, the IMC DAQ system records the displacement and its
frequency from peak-to-peak position to the initial position, and
stores in the computer. The heat sink-piezofan assembly is per-
fectly isolated from external disturbances by means of a large
transparent housing fixed on a rigid platform. For the best perfor-
mance, the piezofans are excited at the resonance frequency
(111 Hz) and 16.5 V. The maximum deflection of each piezofan is
5.01 mm which is symmetric in both directions.
In this study, Central Composite Design (CCD), the most widely
used approach of Response Surface Methodology (RSM), is em-
ployed to optimize d and a with the objective of maximizing tem-
perature drop (DT); the entire DOE procedure is facilitated by the
‘Design Expert 6’ software, similar to the published work [20].

3. Modeling

3.1. Computational domain

The 3D model used in this simulation consisted of a heater, sim-

ple heat sink and three piezofans. Fig. 5 shows the multi-piezofan
mounted in vertical (Fig. 5(a)) and horizontal (Fig. 5(b)) positions.
Fig. 2. Schematic of the heat sink and thermocouple locations. Red area represents
the heat flux imposed (all dimensions in mm and not to scale). (For interpretation of The heat flux employed for numerical simulation was same as that
the references to colour in this figure legend, the reader is referred to the web used in the experiment (9166.7 W/m2).
version of this article.) The size of the computational domain is 120 mm  100 mm 
46.1 mm, which includes heater, heat sink and multi-fans; the
clamp of the piezofan is neglected. The fan boundary is modeled
of the piezofan (PiezoSystem Inc, USA). The length of the un-patch as a moving adiabatic wall whose location in time is set by a
piezofan (lu) is 23 mm. The oscillating drive signal to match the UDF in FLUENT. The fan is modeled as infinitesimally thin wall
resonant frequency of the beam was provided by means of an in- with no thermal conduction through it. No-slip boundary condi-
verter circuit. In the rest of this article a normalized tip gap (d) will tions are applied to the walls such as platform and heat sink. Other
 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525 5517

Fig. 4. Experimental setup to study piezofan deflection.

Fig. 5. Multi-piezofan orientations used in simulation. Fig. 6. The meshed models.

boundaries are treated as pressure boundaries, permitting the air

flow in inward or outward direction. A combination of tetrahedral
and quad literal elements (about half million) is used for meshing
(see Fig. 6).
The fluid domain is divided into two regions namely region ‘a’
and region ‘b’ as shown in Fig. 7. The boundary of these two regions
was imposed as a fluid interface. The purpose of fluid interface is to
avoid damaged fluid cells since region ‘a’ of the fluid domain in-
cludes the deforming beam (piezoelectric fan). The fluid cells can
be greatly skewed due to its movement. This technique can over-
come the problem by flagging the cells from one side of the interface Fig. 7. Mesh generation regions.
5518 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525

as deforming cells. In addition, only single cell spacing is included Differentiating this equation with respect to time gives the velocity
from the tip of the beam to the fluid interface to avoid stretching, of the beam,
which may lead to failure of the fluid cells. Region ‘b’ is distinguished
 
by stationary fluid cells that interact with region ‘a’ through the fluid ðsinðblu Þ  sinhðblu ÞÞðsinðbxÞ  sinhðbxÞÞ
_
wðx; tÞ ¼ Af   xb cosðxb tÞ
interface. The ‘smoothing’ and ‘remeshing’ functions of FLUENT [21] þðcosðblu Þ  coshðblu ÞÞðcosðbxÞ  coshðbxÞÞ
are exploited in the present dynamic meshing. ð8Þ
3.2. Mathematical model
where b values can be calculated from the frequency equation as
The flow is assumed incompressible and turbulent. The govern- follows:
ing equations employed in FLUENT, for describing the transient
fluid flow are as follows: cosðblu Þ  coshðblu Þ ¼ 1 ð9Þ
Continuity:
@q @ This must be solved numerically and yields infinity solutions of b.
þ ðqui Þ ¼ 0 ð1Þ The b value corresponding to the first mode shape is:
@t @xi
Momentum (non-accelerating reference frame):
1:875
b¼ ð10Þ
@ @ @P @ sij lu
ðqui Þ þ ðqui uj Þ ¼  þ þ qg i þ F i ð2Þ
@t @xj @xj @xj
where q is the fluid density, P is the pressure in the fluid, sij is the 3.3. Simulation
viscous stress tensor and gi and Fi are the gravitational acceleration
and external body force in the i-direction, respectively. The first-order upwind discretization scheme is used both for
However, FLUENT allows the user to simulate moving and momentum and energy equations, with the SIMPLE scheme for
deforming domains through the use of user defined function pressure-velocity coupling. For the treatment of the turbulent
(UDF). Dynamic meshes can be used to model flow where the flow at the near-wall region (heat sink surface), the y+ value is
shape of the domain is changing with time due to motion of the set as 1.The beam is assumed to vibrate at a frequency
domain boundaries. The integral form of the transport equation of 100 Hz; this round off (compared to the experimental value
for a general scalar (U), on an arbitrary control volume (V), on a of 111 Hz) is made for numerical stability. Time-step size of
moving mesh is written as: 0.0001 s is chosen for all the cases, with 100 time steps per
Z Z Z Z cycle of fan vibration. This value is arrived after three attempts
d
qUdV þ qUð~
u ~
ug ÞdA ¼ CrUdA þ SU dV ð3Þ to check the trend and proximity with the profiles investigated;
dt V dV dV V for instance, the trials are done by increasing the mesh
where ~u is the flow velocity vector and ~
ug is the grid velocity of the elements (0.34  106  0.6  106) and decreasing the time steps
moving meshes. The first and second terms on the left are the time (2.5  104  5.0  105). With reference to the finest mesh size,
derivative and convective terms respectively. The terms on the right the discrepancy in the temperature of the selected mesh is 0.4%.
are the diffusive and the source terms. C is the diffusion coefficient The total duration of the simulation is selected such that the
and SU is the source term of U. oV represents the boundary of the temperature reaches steady state during this period. It took
control volume V and dA is the area movement. 5000 iterations, which corresponds to approximately three days
Energy: of computation time per each case on a Pentium DualCore proces-
sor (each 2.8 GHz) computer with 2.0 GB of memory. Fig. 8 shows
@ @ @2T the temperature variation over time once the piezofan is actuated.
ðqcp TÞ þ ðqui cp TÞ ¼ k 2 ð4Þ
@t @xj @xj It can be observed that the temperature increases rapidly up to
0.1 s and then slightly increases before attaining steady state at
where cp is specific heat of air, and k is thermal conductivity.
0.2 s. Thus, the quasi-steady temperature can be considered from
The shear stress transport (SST) k–x model as expressed in Eqs.
0.2 s onwards. However, in the current set up these readings are
(5) and (6) is used to describe the flow induced by the piezofan,
taken at 0.5 s, as additional precaution.
which has local turbulence.

@ @ @ @kt e k  Y k þ Sk
ðqkt Þ þ ðqkt ui Þ ¼ þ ðCk ÞþG ð5Þ
@t @xi @xj @xj
@ @ @ @x e x  Y x þ Dx þ Sx
ðqxÞ þ ðqxui Þ ¼ þ ðCx ÞþG ð6Þ
@t @xi @xj @xj

where Ge k represents the generation of turbulence kinetic energy

due to mean velocity gradients, Gx denotes the generation of x,
and Ck and Cx are the effective diffusivity of kt and x, respectively.
Yk and Yx are terms for the dissipations of k and x respectively, due
to turbulence. Dx is the cross-diffusion term, and Sk and Sx are
user-defined source terms.
In present study, the mode shape of a piezoelectric fan is
approximated by that of a clamped-free beam [22]. Assuming a
sinusoidal driving, position of the beam is given by:
 
ðsinðblu Þ  sinhðblu ÞÞðsinðbxÞ  sinhðbxÞÞ
wðx; tÞ ¼ Af   sinðxb tÞ
þðcosðblu Þ  coshðblu ÞÞðcosðbxÞ  coshðbxÞÞ
ð7Þ Fig. 8. Temperature variation over time.
 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525 5519

Table 3
Summary of ANOVA of the quadratic model for DT.

Source Sum of Degree of Mean F-value Prob > F

squares freedom square
Model 10.24 5 2.05 20.56 0.0005
a 1.01 1 1.01 10.14 0.0154
d 7.99 1 7.99 80.18 <0.0001
a2 0.31 1 0.31 3.11 0.1212
d2 0.38 1 0.38 3.78 0.0930
ad 0.14 1 0.14 1.37 0.2797
Residual 0.7 7 0.1
Lack of fit 0.36 3 0.12 1.42 0.611
Pure error 0.34 4 0.084
SD = 0.32 R2 = 0.9363 R2adj ¼ 0:8907 Adeq. Precision = 15.099

Fig. 9. Effect of multi-piezofan on the heat sink temperature.

Table 2
Response values for different experimental conditions.

Run no. Factor A: a Factor B: d Response 1: DT (°C)

1 45.00 0.33 16.60
2 45.00 0.33 15.96
3 90.00 0.33 16.29
4 45.00 0.33 15.94
5 45.00 0.5 15.76
6 45.00 0.33 16.07
7 90.00 0.50 15.53
8 0.00 0.17 17.31
9 45.00 0.17 17.79
10 45.00 0.33 15.67
11 90.00 0.17 18.08
12 0.00 0.50 16.24
13 0.00 0.33 15.53
Fig. 10. Normal probability plot of the studentized residual for heat removal.

In order to ensure that the selected model adequately represents

4. Results and discussion
the real system, the diagnostic plots such as normal probability
plots of the studentized residuals, as well as the predicted versus
4.1. Experimental results
actual value plots, are plotted as shown in Figs. 10 and 11 respec-
tively. It can be observed that the data is almost normally distrib-
4.1.1. Effect of multi-piezofan on temperature
uted even though there is some scattering which is acceptable.
The transient temperature distributions at different locations of
The predicted and experimental values of DT are in good agreement
the heat sink, with and without piezofans, are shown in Fig. 9. First
as shown in Fig. 10. The 3D response surface and contour plots of
of all, keeping the fans off, the heater is turned on (time t = 0), and
the quadratic model as shown in Fig. 12, demonstrate the interac-
the temperature is monitored; the steady state is reached at
tive relationship between d and a on DT. The maximum DT ob-
around t = 80 min. At this condition, the piezofans are turned on
served is 18.2337 °C at d = 0.17 and a = 90°, while the minimum
and a new steady state is reached at around 100 min. The multi-
DT (15.5908 °C) is obtained at d = 0.50 and a = 0°. Thus the optimum
piezofan is seen to cause a maximum of 22% reduction in temper-
values of d and a are 0.17 and 90° respectively, which are further
ature relative to natural convection.
verified by exploiting the optimization facility of Design Expert 6.
4.1.2. Results of DOE analysis
A total of 13 runs obtained from the CCD, and the corresponding 4.1.3. Heat transfer coefficient
responses are shown in Table 2. It is observed that, the tempera- The experimental convective heat transfer coefficient (hhs) in
ture drop (DT) due to the presence of multi-piezofan varies be- each case (without and with the fans) is calculated by using an en-
tween 15.53 °C and 18.08 °C. Table 3 shows the analysis of ergy balance on a control volume between the power input and the
variance (ANOVA) of regression parameters of the predicted re- heat removed from the surface of the heat sink by convection, as
sponse surface quadratic model for DT. The value of correlation proposed by Incropera and Dewitt [23]:
coefficient obtained for DT is R2 = 0.9363. q_ h t h
The final regression model, in terms of their coded factors, is ex- hhs ¼ ð12Þ
T sexp  T 1
pressed by the following second-order polynomial equation (Eq.
(11)). where q_ h is the uniform heat generation (W/m3) from the heater, th
is the heater thickness, and Ts and T1 are the temperatures of the
DT ¼ 15:79 þ 0:41a  1:15d þ 0:33a2 þ 0:37d2  0:18ad ð11Þ heat sink surface and the ambient respectively. Accordingly, hhs
5520 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525

Fig. 11. Comparison of actual and predicted DT.

for three cases namely, best (d = 0.17 and a = 90°), mediocre

(d = 0.33 and a = 45°) and worst (d = 0.50 and a = 0°) are compared
with that of natural convection, as shown in Fig. 13. It is apparent
that, compared to the natural convection case, the multi-piezofan
with vertical orientation and minimum tip gap contributes for
88.5% enhancement in heat transfer coefficient while the mediocre
and worst cases produce 68.8% and 65% respectively. The experi-
mental uncertainties for the parameters (temperature and heat
transfer coefficient) are ±1.72 and 1.88 respectively, which are esti-
mated from 5 sets of readings, by the statistical procedure proposed
by Taylor [24].

4.2. Numerical results

4.2.1. Validation of the structural model

In the present study, the piezofan is approximated as a
clamped-free beam, and Eqs. (7)–(10) are used to define the beam
location during vibration. In order to ensure that this mode shape
matches the desired amplitude, these equations are programmed
in MATLAB software. Fig. 14 demonstrates the successful peak-
to-peak deflection of the beam (piezofan) with respect to the phase Fig. 12. (a) Response surface and (b) contour plots for DT as functions of d and a.
angle (product of angular velocity and time). The mode equations
are then employed into the UDF.

4.2.2. Flow visualization

The transient thermal analysis from FLUENT 6.3.2 is performed
to study the flow field and temperature distribution, for the best
and the worst cases. Fig. 15(a) shows the streamlines, velocity vec-
tors and temperature contours for the best case at t = 0.5 s. The
piezofan produces air streams mainly in the longitudinal and
transverse directions (with respect to the piezofan). The first direc-
tion is generated by the impinged jet, while the second is gener-
ated by the normal force due to the vibration of the beam [25].
So there is fresh cooling air being entrained into the air flow while
the vortices are being formed. As the piezofan vibrates, two regions
of circulation are formed on both sides, forcing the fluid to displace
more rapidly. High velocity regions are found at the front of the
moving beam, which then strikes the heat sink surface and flows
in the outward direction. Fig. 15(b) shows the corresponding
images for the worst case. Unlike the previous case, the velocity
distribution is asymmetric since the states of air on either sides
of the vibrating piezofan are different to each other. With respect Fig. 13. Enhancement in heat transfer coefficient for three different cases.
 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525 5521

Fig. 14. Piezofan deflection with respect to phase angle.

Fig. 15. (a) Predicted streamlines, velocity vectors and temperatures for best case at t = 0.5 s. (b) Predicted streamlines, velocity vectors and temperatures for worst case at
t = 0.5 s.

to the piezofan orientation, the air below the beam is more con- the heater is situated at the bottom center of the heat sink base,
fined and hotter compared to the other side. the central portion of the heat sink gets hotter, and thus the vorti-
Fig. 16 shows the top views of streamlines and vortices at ces generated and concentrated centrally in the first case enables
t = 0.5 s. In the best case (Fig. 16(a)), three vortices are observed ex- efficient heat removal, compared to the second case. Moreover, it
actly at the center of the heat sink surface, while the other case is observed that the magnitude of vorticity for the best case is
(Fig. 16(b)) shows five vortices located away from the center. As 200% greater than that of the other case.
5522 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525

Fig. 17. (a) Temperature contours and (b) labels for best case at t = 0.5 s.

of d and a on the cooling performance. The temperatures at the

leading edge (lhs/lh = 2.76) are almost the same because at this loca-
Fig. 16. (a) Streamlines and vorticity vectors for best case at t = 0.5 s. (b) tion, the induced flows have similar influences for both the cases.
Streamlines and vorticity vectors for worst case at t = 0.5 s. However, temperatures at the center and at the trailing edge are
tremendously different (the best case shows lower temperatures);
this behavior is expected according to the respective flow and vor-
4.2.3. Temperature contours tex profiles.
The temperature contours of the heat sink base, for the best and
worst cases are shown in Figs. 17 and 18 respectively. For the best
case, a rounded rectangular pattern is observed at the heat sink 4.3. Heat transfer coefficient
center (due to the middle fan) and this behavior transforms into
a lobed pattern towards the locations of side fans. A symmetric Fig. 22 shows the heat transfer coefficient (hhs) for the best and
growth of the temperature pattern can be observed in all direc- worst cases. In each case, the heat transfer coefficient decreases as
tions. However, a distinctly different distribution of temperature it moves to the heat sink’s center which is the hottest zone. The in-
is observed in the worst case where the temperature distribution crease in hhs of the best case compared to the worst case is around
is non-uniform. These observations are consistent with the profiles 1% at the leading edge, increases dramatically to 16% at the centre
shown in Figs. 15 and 16 respectively. and then decreases to 10% at the trailing edge. This trend demon-
Figs. 19 and 20 show the corresponding temperature contours strates the crucial roles of a and d on the heat removal
on the fins, at t = 0.5 s. For the first case, a droplet cooling zone is performance.
observed at the centre of the fins and these droplets appear to be
symmetric in the direction of piezofan vibration. The droplet zones 4.4. Comparison of simulation and experimental results
widen laterally due to the induced swirling flow on both sides of
the fan. However, for the other case, asymmetric cooling zone is The predicted and measured temperatures at the locations
observed at one end of the heat sink surface, as the induced flow shown in Fig. 2(b) are compared for the best and worst cases as
freely goes outward. Both profiles agree that this droplet cooling summarized in Tables 4 and 5 respectively. It can be observed that
zone is prominent at the central fins compared to the side fins; this there is good agreement between simulation and experiment, with
is attributed to the enhanced cooling of the central fins due to the a maximum discrepancy of only 1%. Similarly, Fig. 23 shows the
presence of piezofans on both sides. comparison of convective heat transfer coefficients at the location
lhs
Fig. 21 illustrates the comparison of temperature distributions lh
¼ 0 of the heat sink. A maximum of 5% under-prediction is ob-
for the best and worst cases, along the centre line (parallel to the served in the simulation, which is attributed to the omission of
x-direction referring to Fig. 2) of the heat sink surface. The signifi- radiation effect in the numerical model; this presumption is
cant difference between the curves indicates the critical influence confirmed by deducting the experimentally determined radiative
 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525 5523

Fig. 20. Temperature contours of the fins for the worst case at t = 0.5 s.

Fig. 18. (a) Temperature contours and (b) labels for worst case at t = 0.5 s.

Fig. 21. Temperature plots on the reference line for the comparison of best and
worst cases.

Fig. 19. Temperature contours of the fins for the best case at t = 0.5 s.

contribution, from the predicted value. Also, the thickness of the

thermal compound in the simulation is treated as uniform in all
the cases; the real thickness was not necessarily so, as the com-
pound was manually applied. Furthermore, the ambient tempera-
ture in the experiments is measured with a single thermocouple,
whereas the average fluid temperature inside the domain is used
in the simulation. In addition, the lack of match between the exper-
imental and predicted frequencies might contribute to the incon- Fig. 22. Heat transfer coefficient plots on the reference line for the comparison of
sistency which leads to difference in performance. best and worst cases.
5524 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525

Table 4
Comparison of experimental and simulation temperatures (best case).

Temperature points Ts 1 2 3 4 5 6 7
Exp. (K) 317.48 312.34 310.64 312.06 306.46 306.98 306.83 306.33
a = 90 (K) 317.80 312.82 311.66 312.44 310.76 311.55 311.29 309.39
DTexp–sim (K) 0.32 0.48 1.02 0.38 4.3 4.57 4.46 3.06
%DTexp–sim 0.1 0.15 0.33 0.12 1.4 1.49 1.45 1.0
% Average DTexp–sim 0.76

Table 5
Comparison of experimental and simulation temperatures (worst case).

Temperature points Ts 1 2 3 4 5 6 7
Exp. (K) 320.39 312.76 312.82 311.89 306.08 310.64 305.31 308.06
a=0 (K) 321.52 314.44 313.94 314.02 311.94 313.01 310.79 313.09
DTexp–sim (K) 1.13 1.68 1.12 2.13 5.86 2.37 5.48 5.03
%DTexp–sim 0.35 0.54 0.36 0.68 1.91 0.76 1.79 1.63
% Average DTexp–sim 1.0

W. Amri for their contributions in fabrication of the test apparatus

and in conducting the noise and vibration experiments. Special
thanks also dedicated to M. Hashim for helping in conducting the
control and instrumentation experiments.

References

[1] T. Acikalin, S.V. Garimella, Analysis and prediction of the thermal performance
of piezoelectrically actuated fans, Heat Transfer Eng. 30 (6) (2009) 487–498.
[2] M. Toda, Theory of air flow generation by a resonant type PVF2 bimorph
cantilever vibrator, Ferroelectrics 22 (1979) 911–918.
[3] M. Toda, Voltage-induced large amplitude bending device-PVF2 Bimorph—its
properties and applications, Ferroelectrics 32 (1981) 127–133.
[4] A. Ihara, H. Watanabe, On the flow around flexible plates, oscillating with large
amplitude, J. Fluids Struct. 8 (1994) 601–619.
[5] T. Acıkalın, A. Raman, S.V. Garimella, Two-dimensional streaming flows
induced by resonating thin beams, J. Acoust. Soc. Am. 114 (No. 4-Pt. 1)
(2003) 1785–1795.
[6] M.K. Abdullah, M.Z. Abdullah, M.V. Ramana, C.Y. Khor, K.A. Ahmad, M.A.
Mujeebu, Y. Ooi, Z. Mohd Ripin, Numerical and experimental investigations on
effect of fan height on the performance of piezoelectric fan in microelectronic
Fig. 23. Comparison of predicted and experimental heat transfer coefficients for cooling, Int. Commun. Heat Mass Transfer 36 (2009) 51–58.
best and worst cases. [7] M.K. Abdullah, M.Z. Abdullah, M.V. Ramana, C.Y. Khor, K.A. Ahmad, M.A.
Mujeebu, Y. Ooi, Z. Mohd Ripin, Effect of piezoelectric fan height on flow and
heat transfer for electronics cooling applications, in: Proceedings of the EMAP
2008 Conference, Taipei, Taiwan, 2008, pp. 165–170.
5. Conclusions [8] R.R. Schmidt, Local and average transfer coefficients on a vertical surface due
to convection from a piezoelectric fan, in: Proceedings of the I-THERM IV
Optimization of tip gap and orientation of piezofans is per- Conference, Washington, 1994, pp. 41–49.
[9] T. Acıkalın, S.M. Wait, S.V. Garimella, A. Raman, Experimental investigation of
formed by using the DOE approach, to maximize the heat removal
the thermal performance of piezoelectric fans, Heat Transfer Eng. 25 (1) (2004)
performance of finned heat sink. A 3D simulation using FLUENT 4–14.
6.3.2 is carried out to better understand the induced flow and [10] T. Acıkalın, I. Sauciuc, S.V. Garimella, Piezoelectric actuators for low-form
factor electronics cooling, in: Proceedings of the IPACK2005-ASME InterPACK
the resulting heat transfer. The turbulence was incorporated by
Conference, San Fransisco, CA, 2005, pp. 939–943.
the shear stress transport (SST) k–x model. The experimental re- [11] J.H. Yoo, J.I. Hong, W. Cao, Piezoelectric ceramic bimorph coupled to thin metal
sults are well matched with the predictions. By the use of DOE, sig- plate as cooling fan for electronic devices, Sensors Actuat. 79 (2000) 8–12.
nificant savings in effort and time are achieved to conduct the [12] H. Hosaka, K. Itao, Coupled vibration of microcantilever array induced by
airflow force, ASME J. Vib. Acoust. 124 (2002) 26–32.
required experiments and analyze the coupled influence of tip [13] S. Basak, A. Raman, Hydrodynamic coupling between micromechanical beams
gap and piezofan orientation. It is observed that multi-piezofan oscillating in viscous fluids, Phys. Fluids 19 (1017105) (2007) 1–13.
with vertical orientation and a tip gap d = 0.17 yields 88% enhance- [14] M. Kimber, R. Lonergan, S.V. Garimella, Experimental study of aerodynamic
damping in arrays of vibrating cantilevers, J. Fluids Struct. 25 (2009) 1334–
ment in the convective heat transfer coefficient compared to natu- 1347.
ral convection. The present finding would be a valuable input in [15] M. Kimber, S.V. Garimella, Cooling performance of arrays of vibrating
the application of finned heat sink with multi-peizofan in micro- cantilevers, J. Heat Transfer 131 (2009) 1–8.
[16] M.K. Abdullah, N.C. Ismail, M.Z. Abdullah, K.A. Ahmad, M.A. Mujeebu, F. Husin,
electronic cooling. This study may be extended for lesser tip gaps, Z. Mohd Ripin, Heat transfer enhancement using piezoelectric fan in electronic
and with various piezofan-heat sink configurations. cooling-experimental and numerical observations, ISI Bilimi ve Teknigi
Dergisi-journal of thermal science and technology 32 (2012) 1.
[17] T. Acikalin, S.V. Garimella, A. Raman, J. Petroski, Characterization and
optimization of the thermal performance of miniature piezoelectric fans, Int.
Acknowledgements J. Heat Fluid Flow 28 (2007) 806–820.
[18] J. Petroski, M. Arik, M. Gursoy, Optimization of piezoelectric oscillating fan-
This work is supported by the Fundamental Research Grant cooled heat sinks for electronics cooling, IEEE Trans. Comp. Pack. Technol. 33
(1) (2010) 25–31.
Scheme (FRGS) under contract 1001/PMEKANIK/811076. The [19] M.K. Abdullah, N.C. Ismail, M.Z. Abdullah, M. Abdul Mujeebu, K.A. Ahmad, Z.
authors would like to acknowledge S.F. Shaker, T.Y. Chong, and Mohd Ripin, Effects of tip gap and amplitude of piezoelectric fans on the
 M.K. Abdullah et al. / International Journal of Heat and Mass Transfer 55 (2012) 5514–5525 5525

performance of heat sinks in microelectronic cooling, Heat Mass Transfer [23] F.P. Incropera, D.P. Dewitt, Fundamentals of Heat and Mass Transfer, fourth
(2012), http://dx.doi.org/10.1007/s00231-011-0944-z. ed., John Wiley & Sons, Singapore, 1996. pp. 103–105.
[20] H.A. Mohamed, M.Z. Elnaggar, M. Abdullah, Abdul Mujeebu, Experimental [24] J.R. Taylor, An Introduction to Error Analysis: The Study of Uncertainties in
investigation and optimization of heat input and coolant velocity of finned Physical Measurements, second ed., vol. 94965, University Science Books,
twin U-shaped heat pipe for CPU cooling. Experimental Techniques, http:// Sausalito, CA, 1997.
dx.doi.org/10.1111/j.1747-1567.2011.00757.x. [25] S.F. Liu, R.T. Huang, W.J. Sheu, C.C. Wang, Heat transfer by a piezoelectric fan
[21] Modeling Flows Using Sliding and Deforming Meshes, FLUENT 6.3.26, FLUENT on a flat surface subject to the influence of horizontal/vertical arrangement,
Inc. Lebanon, New Hampshire, 2006, Chapter 11.3. International Journal of Heat and Mass Transfer 52 (2009) 2565–2570.
[22] L. Meirovitch, Principle and Techniques of Vibration, Prentice Hall, Englewood
Cliffs, Singapore, 1997.