Considerations on Optimum Design of Micro Heat Pipe Sinks Using Water as
Working Fluid
Except where reference is made to the work of others, the work described in this
dissertation is my own or was done in collaboration with my advisory committee. This
dissertation does not include proprietary or classified information.
Florentina Simionescu
Certificate of Approval:
Roy W. Knight
Assistant Professor
Mechanical Engineering
Daniel K. Harris, Chair
Associate Professor
Mechanical Engineering
Jay M. Khodadadi
Professor
Mechanical Engineering
Amnon J. Meir
Professor
Mathematics and Statistics
Joe F. Pittman
Interim Dean
Graduate School
Considerations on Optimum Design of Micro Heat Pipe Sinks Using Water as
Working Fluid
Florentina Simionescu
A Dissertation
Submitted to
the Graduate Faculty of
Auburn University
in Partial Fulfillment of the
Requirements for the
Degree of
Doctor of Philosophy
Auburn, Alabama
December 15, 2006
Considerations on Optimum Design of Micro Heat Pipe Sinks Using Water as
Working Fluid
Florentina Simionescu
Permission is granted to Auburn University to make copies of this dissertation at its
discretion, upon the request of individuals or institutions and at
their expense. The author reserves all publication rights.
Signature of Author
Date of Graduation
iii
Vita
Florentina Simionescu, daughter of Gheorghe and Ioana Popa, was born on March
22, 1973, in Br?aila, Romania. She graduated from Bogdan Petriceicu Ha?sdeu High School
in Buz?au, in 1991. In the same year she entered Transilvania University of Bra?sov and
graduated in 1996 with a Bachelor of Science degree in Mechanical Engineering, major-
ing in Mechatronics. She received a Master of Science degree in Mechanical Engineering
and a Master of Environmental Engineering in 1997 and 1999, respectively from the same
university. Between 1998-2000 she worked as an instructor in the Department of Preci-
sion Mechanics, Transilvania University. In August 2000 she joined the graduate program
of the Mathematics Department at Auburn University and received a Master of Applied
Mathematics degree in May 2005. In January 2001 she started her Ph.D. studies in the
Department of Mechanical Engineering at Auburn University. She married Petru Aurelian
Simionescu on February 26, 2000.
iv
Dissertation Abstract
Considerations on Optimum Design of Micro Heat Pipe Sinks Using Water as
Working Fluid
Florentina Simionescu
Doctor of Philosophy, December 15, 2006
(M.A.M., Auburn University, 2005)
(M.Env.Eng., Transilvania University of Bra?sov, 1999)
(M.S., Transilvania University of Bra?sov, 1997)
(B.S., Transilvania University of Bra?sov, 1996)
133 Typed Pages
Directed by Daniel K. Harris
The successful operation of RF wideband gap devices dissipating high levels of heat
flux requires effective cooling techniques. One of the most promising thermal management
strategies is the use of micro heat pipes (MHP). These devices are very thin profile heat
spreaders that can be directly attached to the GaAs or GaN substrate whose function is
to allow the spreading of the heat flux almost laterally within a 300 ?m thickness. The
objective of the work presented was to estimate the convective heat transfer coefficient
of a micro-channel heat sink corresponding to a maximum amount of heat removed from
heat source placed on the top surface of the sink. This approach taken used an optimal
control technique in which the solution of the heat equation is controlled by the convective
boundary condition by taking the heat transfer coefficient as the control parameter. A
conjugate gradient method was used to solve the optimal control problem. The results
show that the temperature distributions corresponding to the controlled solution are lower
than those correspondingto the uncontrolled solution. Thedifference between the controlled
v
and uncontrolled temperatures (4K?10K) is smaller than in the case of a planar spreader.
This suggests that the convective coefficient corresponding to a finned spreader approaches
by design an optimum distribution.
The effect of liquid charge on the performance of MHP sinks is important. A too
large amount of fluid leads to a condenser flooding, while a too low fluid charge leads to
an evaporator dry-out and an increase in channel wall temperature. An iterative scheme
was devised to compute the liquid charge corresponding to the maximum heat transport
capacity of the pipe.
This study can provide guidance in designing MHP sinks, which have emerged as an
effective technique for cooling electronic components.
vi
Acknowledgments
I would like to express my deepest gratitude to my advisor, Dr. Daniel Harris, Associate
Professor of Mechanical Engineering, for his guidance and support toward the completion
of this dissertation.
Special thanks are extended to Dr. A.J. Meir, Professor of Mathematics, for his invalu-
able help when needed mostly. I would also like to acknowledge Dr. Roy Knight, Assistant
Professor of Mechanical Engineering, and Dr. Jay Khodadadi, Professor of Mechanical
Engineering for their help and suggestions as committee members, as well as Dr. William
Ashurst, the outside reader, for his thorough review of the manuscript.
I am grateful to my parents and sister for their love and moral support. Last but not
least, I would like to thank my husband, Aurelian for his love, patience, and permanent
encouragement that made this work possible.
In loving memory of my grandmother Caterina
vii
Style manual or journal used International Journal of Mass and Heat Transfer(together
with the style known as ?aums?).
Computer software used The document preparation package TEX (specifically LATEX)
together with the departmental style-file aums.sty.
viii
Table of Contents
List of Figures xi
1 Introduction 6
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2 Optimal convective heat transfer coefficient of a planar spreader 24
2.1 Problem specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.2 Existence and uniqueness of solution . . . . . . . . . . . . . . . . . . . . . . 28
2.3 Optimal control for model problem . . . . . . . . . . . . . . . . . . . . . . . 29
2.3.1 The direct problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.3.2 The optimal control problem . . . . . . . . . . . . . . . . . . . . . . 31
2.3.3 Existence and uniqueness of an optimal control . . . . . . . . . . . . 32
2.3.4 Conjugate gradient algorithm . . . . . . . . . . . . . . . . . . . . . . 33
2.3.5 Sensitivity problem and calculation of the search step size . . . . . . 34
2.3.6 The adjoint problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.3.7 Computational algorithm . . . . . . . . . . . . . . . . . . . . . . . . 38
2.4 Implementation and numerical procedure . . . . . . . . . . . . . . . . . . . 39
2.5 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3 Optimal design of a micro heat pipe 44
3.1 Problem specification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2 The direct problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2.1 Heat transfer equations corresponding to the GaAs slab and silicon
wafer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2.2 1-D Micro heat pipe model . . . . . . . . . . . . . . . . . . . . . . . 49
3.3 The optimal control problem with Bi and geometrical parameters as controls 62
3.3.1 Conjugate gradient algorithm . . . . . . . . . . . . . . . . . . . . . . 63
3.3.2 Sensitivity problem and calculation of the search step sizes . . . . . 64
3.3.3 Adjoint problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.3.4 Computational algorithm . . . . . . . . . . . . . . . . . . . . . . . . 75
3.3.5 Implementation and numerical procedure . . . . . . . . . . . . . . . 76
3.3.6 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . 76
3.4 Liquid charge corresponding to maximum heat transport capacity of the pipe 78
3.4.1 Computational procedure . . . . . . . . . . . . . . . . . . . . . . . . 79
3.4.2 Results and discussions . . . . . . . . . . . . . . . . . . . . . . . . . 80
ix
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4 Appendix 86
x
List of Figures
2.1 Computational domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.2 Structure of matrix A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.1 Schematic of a MHP sink. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2 Computational domain. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.3 Schematic of a MHP. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4 Control volume used to derive the conservation equations . . . . . . . . . . 52
3.5 A cross section of the MHP with liquid recessed into the corner . . . . . . . 53
3.6 Conservation of mass . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
3.7 Conservation of momentum . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1 Graphs of controlled and uncontrolled solutions for different values of ?. . . 87
4.2 Distributions of Bi for different values of ?. . . . . . . . . . . . . . . . . . . 88
4.3 Distributions of T ?T? at the chip-spreader interface for different values of ?. 89
4.4 Distributions of T?T? through a vertical symmetry plane for different values
of ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
4.5 Distributions of optimal Bi corresponding to N = 60 and different channel
depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.6 Distributions of optimal Bi corresponding to N = 58 and different channel
depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.7 Distributions of optimal Bi corresponding to N = 56 and different channel
depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
4.8 Distributions of optimal Bi corresponding to N = 30 and different channel
depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
xi
4.9 Temperature distributions T ?T? corresponding to N = 60. . . . . . . . . 95
4.10 Temperature distributions T ?T? corresponding to N = 58. . . . . . . . . 96
4.11 Temperature distributions T ?T? corresponding to N = 56. . . . . . . . . 97
4.12 Temperature distributions T ?T? corresponding to N = 30. . . . . . . . . 98
4.13 Distribution of the meniscus radius for different values of ? and channel
depths (N = 60). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
4.14 Distributions of the liquid and vapor velocities for H1 = 75?m and different
values of ? (N = 60). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
4.15 Distributions of the liquid and vapor velocities for H1 = 100?m and different
values of ? (N = 60). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
4.16 Distributions of the liquid and vapor velocities for H1 = 125?m and different
values of ? (N = 60). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
4.17 Distributions of the liquid and vapor velocities for H1 = 150?m and different
values of ? (N = 60). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
4.18 Distributions of the liquid and vapor pressures for different values of ? and
channel depths (N = 60). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
4.19 Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding
to the maximum heat transport capacity of the pipe, for N = 60 and different
channel depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.20 Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding
to the maximum heat transport capacity of the pipe, for N = 58 and different
channel depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
4.21 Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding
to the maximum heat transport capacity of the pipe, for N = 56 and different
channel depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
4.22 Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding
to the maximum heat transport capacity of the pipe, for N = 30 and different
channel depths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
4.23 Variations of the liquid charge with r0 for different channel depths (N = 60). 109
xii
4.24 Variations of MHP heat transport capacity with r0 for different channel
depths (N = 60). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
4.25 Variations of liquid charge with channel depth for different values of ?. . . . 111
4.26 Percent of volume charge corresponding to maximum heat transport capacity
of the pipe. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
4.27 Variations of MHP heat transport capacity with channel depth for different
values of ?. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
xiii
Nomenclature
a bilinear form
cl,clw,ci geometrical constants (functions of ?, ?, and r)
dc height of chip
fl friction coefficient for liquid phase
fv friction coefficient for vapor phase
h convective heat transfer coefficient
hfg latent heat of vaporization
kc coefficient of thermal conductivity for chip
ks coefficient of thermal conductivity for planar spreader
kSi coefficient of thermal conductivity for Si
kGaAs coefficient of thermal conductivity for GaAs
l width of channel
lc length of chip
p channel aspect ratio
qr maximum heat transport capacity of MHP as a fraction of the input heat
qv volumetric heat source
qw heat flux entering the MHP
r meniscus radius
rmax maximum radius of curvature
w1 length of heating area
w2 length of GaAs slab
wc width of chip
1
Ai liquid-vapor interfacial area inside MHP
Al liquid cross sectional area inside MHP
Av vapor cross sectional area inside MHP
Alw wall/liquid area inside MHP
Avw wall/vapor area inside MHP
Bi Biot number
Bi? minimum Biot number
Bo Bond number
C0,1 space of Lipschitz continuous functions
Ca capillary number
Ds height of planar spreader
DH hydraulic diameter of empty pipe
DHl hydraulic diameter of liquid phase
DHv hydraulic diameter of vapor phase
H1 (?) Sobolev space
H1 height of fin (channel)
H2 height of Si substrate
H3 height of GaAs slab
J functional to be minimized
L width of computational domain
L2 (?) space of square integrable functions
La length of MHP adiabatic section
Le length of MHP evaporator section
2
Ls length of planar spreader
Lt width of MHP heat sink
N number of channels
P,P1,P2 search directions in conjugate gradient algorithm
Pl liquid pressure
Pv vapor pressure
Pref reference pressure
?Q dissipated power
RT total thermal resistance of a heat sink
Rel Reynolds number for liquid phase
Rev Reynolds number for vapor phase
S chip/planar spreader interfacial area
T? ambient temperature
T? minimum temperature
Tc temperature distribution within the chip
Tf temperature of coolant
Ts temperature distribution within the spreader
Tw micro-channel wall temperature
TA temperature distribution within GaAs slab
TB temperature distribution within Si base plate
TC temperature distribution within separating wall (fin)
Uad set of admissible controls
Ul liquid velocity
3
Uv vapor velocity
Uref reference velocity
Vil liquid phase interfacial velocity
Viv vapor phase interfacial velocity
Ws width of planar spreader
W length of MHP sink
Greek letters
? weighting coefficient
?,?1,?2 search sizes in conjugate gradient algorithm
?,?1,?2 conjugate coefficient in conjugate gradient algorithm
?? boundary of domain ?
?1,?2 adjoint variables corresponding to chip and spreader respectively
? wall/liquid contact angle
?1 MHP sink tilt angle
?1,?2,?3 adjoint variables associated with p
?1,?2,?3 adjoint variables associated with Bi
?l liquid density
?v vapor density
? surface tension
?li interfacial shear stress for liquid phase
?vi interfacial shear stress for vapor phase
?lw liquid/wall shear stress
4
?vw vapor/wall shear stress
? pipe corner angle
? micro-channel cross section
?s bottom surface of planar spreader
?p small perturbation of p
?Bi small perturbation of Bi
?Tc sensitivity temperature within the chip
?Ts sensitivity temperature within the planar spreader
?T1A sensitivity temperature within GaAs slab when p is perturbed by ?p
?T1B sensitivity temperature within Si base plate when p is perturbed by ?p
?T1C sensitivity temperature within fin when p is perturbed by ?p
?T2A sensitivity temperature within GaAs slab when Bi is perturbed by ?Bi
?T2B sensitivity temperature within Si base plate when Bi is perturbed by ?Bi
?T2C sensitivity temperature within fin when Bi is perturbed by ?Bi
?TS?C temperature rise above the input coolant temperature
? computational domain for planar spreader
?c chip subdomain
?s planar spreader subdomain
5
Chapter 1
Introduction
1.1 Background
With the increase in heat dissipation from microelectronic devices and their size reduc-
tion, thermal management becomes a more important element of electronic product design.
Both the reliability and life expectancy of electronic equipment are inversely related to the
component temperature of the equipment. The relationship between the reliability and the
operating temperature of a typical silicon semi-conductor device shows that a reduction in
the temperature corresponds to an exponential increase in the reliability and life expectancy
of the device [35], [10], [75], [56]. Therefore, long life and reliable performance of a com-
ponent may be achieved by effectively controlling the device operating temperature within
the limits set by the device design engineers.
Heat sinks are devices that enhance heat dissipation from a hot surface, usually the case
of a heat generating component, to a cooler ambient. The primary purpose of a heat sink
is to maintain the device temperature below the maximum allowable temperature specified
by the device manufacturers.
When designing or selecting an appropriate heat sink one needs to examine various
parameters that affect not only the heat sink performance itself, but also the overall perfor-
mance of the system. The performance of a heat sink is measured by its thermal resistance
RT = ?TS?C/ ?Q (1.1)
6
where ?TS?C is the temperature rise of the electronic component above the input coolant
temperature and ?Q is the dissipated power. The total thermal resistance is a sum of three
components that account for conduction through the silicon (aluminum/copper) substrate,
convection from the substrate to the cooling fluid, and resistance due to the heating of fluid
as it absorbs the energy passing through the substrate, respectively. The optimization of
a heat sink aims to minimize the total thermal resistance. The optimization procedure is
performed subject to design constraints such as:
? Available pressure drop;
? Cross sectional geometry of incoming flow;
? Ambient fluid temperature;
? Amount of required heat dissipation;
? Maximum size of the heat sink;
? Orientation with respect to gravity, etc.
Usually the sinks are equipped with extended surfaces (fins) that increase the amount of
heat removed. The parameters over which the designer has control when performing the
optimization of a heat sink typically include:
? Fin height;
? Fin length;
? Fin thickness;
? Number/density of fins;
7
? Fin shape/profile;
? Base plate thickness;
? Heat sink material;
? Cooling fluid properties, etc.
There is always an interdependence among the above mentioned parameters. The impact
one parameter has on the performance of a heat sink can not be determined without con-
currently considering the contribution of the other parameters. For example, a greater fin
height provides additional surface area for heat dissipation and improves the overall thermal
performance. However, if the available volumetric flow rate is fixed, the overall performance
may deteriorate by increasing the fin height.
1.2 Objectives
The purpose of this work is to study the performance of a heat sink equipped with
an array of micro heat pipes. The geometrical parameters of the heat sink and the heat
transfer coefficient are optimized using an optimal control technique. These optimal design
parameters correspond to a maximum amount of heat removed from a heat source placed
on top of the heat sink. The liquid charge corresponding to the maximum heat transport
capacity of the pipe is then computed.
1.3 Literature review
During the last two decades, several cooling schemes ranging from simple planar ther-
mal spreaders [9], [19], to more sophisticated two-phase flow heat sinks and microjet cooling
8
devices have been proposed and extensively studied both experimentally and theoretically.
Some of the most representative works are briefly presented in the sequel.
Peles et al. [48] analyzed the performance of micro pin fin heat sinks. A simplified
expression for the total thermal resistance was derived and experimentally validated. Geo-
metrical and thermo-hydraulic parameters affecting the total thermal resistance have been
discussed. It has been found that very low thermal resistances are achievable using a pin fin
heat sink. The thermal resistances achieved using micro pin fin heat sinks were as low as the
ones with micro-channel convective flows. The advantage of using pin fin configurations was
the design flexibility in the geometrical selection of the pin shapes and their spacing. Exper-
iments were performed by Wei and Honda [71] to study the effects of square micro pin fins
on boiling heat transfer from silicon chips immersed in a pool of degassed or gas-dissolved
FC-72. The micro pin finned chips showed a considerable heat transfer enhancement in the
nucleate boiling region and increase in the critical heat flux, as compared to the smooth
chip. The maximum value of allowable heat flux (84.5W/cm2), 4.2 times as large as that for
the smooth chip, was obtained for a fin height of 270?m, a fin width of 50?m, and a liquid
subcooling of 45K. Another numerical study on pin fin heat exchangers was carried out by
Saha and Acharya [55] to analyze the unsteady 3-D flow and heat transfer in a parallel-plate
channel heat exchanger with in line arrays of periodically mounted rectangular cylinders
(pins) at various Reynolds numbers and geometrical configurations. Results were presented
to highlight different parametric effects, using both the instantaneous snapshots of the flow
and heat transfer fields and the averaged (space- and time-averaged in case of unsteady
flows) integral parameters, such as the friction factor and the Nusselt number. It was found
that the thermal performance increases significantly when the flow becomes unsteady.
9
The thermal performance of liquid cooling heat sinks is usually one or two orders of
magnitude greater than the one of the heat sinks using air for cooling. The air cooling
devices on the other hand, are cheap, easy to implement, and reliable. Micro-channel heat
exchangers combine the attributes of very high surface area to volume ratio, large convective
heat transfer coefficient, small mass and volume, and small coolant inventory. Numerous
investigations have been conducted on forced convection in micro-channels concerning both
single phase and two-phase flows. A comparison of single-phase and two-phase (liquid-
vapor) heat sinks shows that there is a significant enhancement in heat transfer in the
latter case. The latent heat transfer in the two-phase systems allows the same amount
of heat transfer at a much lower temperature difference (between the evaporator and the
condenser).
The micro-channel heat sink concept was first introduced by Tuckerman and Pease
in the early 1980s [68]. They performed an optimization of the micro-channel heat sink
design subject to a fully developed and laminar in nature flow through the channels and
a laminar Nusselt number. The channel to fin width ratio, the pressure drop through
the fin array, the pumping work, the planar dimensions and the fin efficiency were fixed.
They fabricated and tested several micro-channel heat sinks in a 1cm?1cm silicon wafer.
The experimental findings were in good agreement with the theoretical predictions. The
heat sink with channels separated by 50?m thick walls, having a width of 50?m and a
depth of 302?m proved to be the optimal design solution. Using water as cooling fluid, the
micro-channel heat sink was capable of dissipating 790W/cm2 with a maximum substrate
temperature raise of 70?C above the water inlet temperature for a pressure drop of 31psi.
10
Following the work of Tuckerman and Pease [68] much research has been conducted for
micro-channel heat sinks.
Weisberg et al. [72] analyzed the micro-channel heat exchangers solving for the coupled
temperature fields within the solid substrate and the working fluid. A 2-D heat transfer was
considered in both the fins and the fluid charged channels. The minimum thermal resistance
was attained for a fin to channel width ratio equal to 1, this result being in accordance with
the one obtained by Tuckerman and Pease [68].
Knight et al. [32] performed a micro-channel heat sink optimization allowing the relax-
ation of the aspect ratio, the channel to fin width ratio, Nusselt number, flow regime and
volumetric flow rate. They also included in their analysis the effects of the developing chan-
nel flow. The results indicated that when the pressure drop through the channel is small,
laminar solutions yield lower thermal resistance than turbulent solutions. Conversely, when
the pressure drop is large the optimal thermal resistance is found in the turbulent region.
With the relaxation of the above constraints, configurations that produced significant im-
provement in thermal resistance over that presented by Tuckerman and Pease [68] were
found. When the turbulent regime was allowed, the thermal resistance was reduced by 35%
from that of Tuckerman and Pease [68]. These theoretical predictions were experimentally
verified [33]. Three heat sink systems from aluminum alloy with 5, 11 and 8 fins, respec-
tively were built and tested, using air as coolant. The fin system with the lowest thermal
resistance was predicted to occur for 9 channels, or 8 fins, and a ratio of fin to channel width
of 0.565. This configuration was predicted to yield a turbulent flow, and a lower resistance
than any laminar configuration operating under the same constraints. At 50W of power
11
dissipation, the optimal design operated about 8?C cooler than the 5 fin system, and about
4?C cooler than the 11 fin system.
The entrance effects on the thermal performance of the micro-channel heat sink were
also examined by Ryu et al. [53]. They performed a 3-D heat transfer analysis which was
further incorporated in an optimization scheme to find the optimal design parameters of
the micro-channel heat sink, that minimize the thermal resistance subject to a specified
pumping power and a channel aspect ratio less than 10. They pointed out that the thermal
entrance effect is substantial when the working fluid is water (Pr ? 7). They varied the
number of channels and obtained the associated design variables and the thermal resistance.
A low sensitivity of the thermal resistance to the number of channels was found. When the
number of channels was reduced to half, the thermal resistance increased by about 15% of
the optimal value. Itwas shown that boththe optimal dimensionsand thethermal resistance
have a power-law dependence on the pumping power. The optimal design parameters for
a pumping power of 2.56W and an input heat flux of 100W/cm2 were found as follows:
the number of channels N = 124, the channel width W = 45.3?m, and the channel height
H = 453?m for a silicon substrate of 100?m. The correspnding thermal resistance was
0.069. The same authors studied the performance of the manifold micro-channel heat sinks
[54]. These heat sinks have an array of manifold dividers perpendicular to the micro-
channels and placed atop of them. The coolant flows through the alternating inlet and
outlet manifolds in the direction normal to the heat-sink base, to and from the segmented
micro-channels, reducing the flow path to a small fraction of the total length of a heat sink.
The shortened flow path reduces the pressure drop and restrains the growth of the thermal
boundary layer in the streamwise direction. Compared to the traditional heat sink for the
12
identical heat load and pumping power, the thermal resistance was lowered by more than
50% while the maximum temperature variation on the heated wall was improved by tenfold.
Peng and Peterson [50] investigated experimentally the single-phase forced convective
heat transfer and flow characteristics of water in micro-channel structures/plates with small
rectangular channels having hydraulic diameters of 0.133mm?0.367mm and different aspect
ratios. The laminar heat transfer was found to be dependent on the aspect ratio H/W and
the ratio of the hydraulic diameter to the center-to-center distance of the micro-channels.
The turbulent heat transfer on the other hand was found to be a function of a new dimen-
sionless variable, Z = min(H,W)/max(H,W), such that Z = 0.5 will give the optimum
configuration for turbulent heat transfer regardless of the groove aspect ratio. The friction
factor or flow resistance reached a minimum value as Z approached 0.5. Empirical corre-
lations were suggested for calculating both the heat transfer and pressure drop. Another
experimental work by Peng and Peterson [49], [51] showed that the range of the transition
zone, and the heat transfer characteristics of both the transition and laminar flow regimes,
were strongly affected by the liquid temperature, liquid velocity and micro-channel size.
Zhao and Lu [76] used two approaches in analyzing the micro-channel heat sinks: the
porous medium model and the fin model. In the porous medium approach, the modified
Darcy equation for the fluid and the two-equation model for the heat transfer between
the solid and fluid phases were employed. It was shown that the overall Nusselt number
increases with increasing aspect ratio and decreases with increasing effective thermal con-
ductivity ratio. Whereas the porous medium model predicted the existence of an optimal
porosity, the fin approach predicted that the heat transfer capability of the heat sink in-
creases monotonically with the porosity (defined as the ratio of the void volume to the total
13
volume). The effect of turbulent heat transfer within the micro-channel was also discussed
and it was found that turbulent heat transfer results in a decrease of optimum porosity in
comparison to that for laminar flow. The concept of micro-channel cooling in combination
with micro heat pipes was proposed. Using a lumped capacitance technique, the micro heat
pipes were modeled as homogeneous, solid regions embedded in the plates with an effective
thermal conductivity that was assumed to be approximately 10 times that of silicon or cop-
per. It was observed that the enhancement in heat transfer due to the incorporation of heat
pipes becomes significant as the channel aspect ratio increases and for highly conducting
coolants, such as water. Another porous medium model developed by Wang et al. [70] for
two phase flows in mini-channels was used to predict the capillary limit of an operating
micro heat pipe. Imke [26] developed a code to simulate the micro-channel flow and heat
transfer in compact heat exchangers. The method was based on a forced convection porous
medium approach combined with conventional pipe flow closure relations and solved for
outlet temperatures, pressure losses, and vapor volume fractions. The code was applied to
single phase cross flow and counter flow heat exchangers using water as working fluid. The
experimental findings were in good agreement with the theoretical predictions.
Vafai and Zhu [69] proposed the design of two-layer micro-channel heat sinks based
on stacking two layers of micro-channels one atop the other, with coolant flow in opposite
directions in each of the micro-channel layers. This solution aimed to reduce the undesirable
temperature gradients in the streamwise direction. The maximum temperature difference
along the pipe was found to be 3 times lower than that of a one-layer micro-channel system.
Some optimization issues for design parameters were addressed. The optimum ratio of
channel width to fin width was found to be 0.6. In the analysis of one-layer structure
14
was pointed out that smaller channel size and higher aspect ratio can reduce the thermal
resistance. This observation proved to hold true for two-layer structures as well. For the
ratio of the pressure drop in the channel array closest to the heat source to the pressure
drop in the second array the recommended optimum value was 2.5?3. Another parameter
discussed was the channel length and was shown that a longer micro-channel design can
indeed substantially reduce the streamwise temperature variation.
Wu and Cheng investigated the effect of surface condition on the performance of silicon
micro-channels [73]. Based on 168 experimental data points, dimensionless correlations for
the Nusselt number and the apparent friction constant were obtained for the flow of water
in trapezoidal micro-channels, having different geometric parameters, surface roughness and
surface hydrophilic properties. It was found that the laminar Nusselt number and apparent
friction constant of the trapezoidal micro-channels increase with the increase of surface
roughness. This increase is more obvious at large Reynolds numbers than at low Reynolds
numbers. The Nusselt number and apparent friction constant of the trapezoidal micro-
channels having strong hydrophilic surfaces (thermal oxide surfaces) were larger than those
having weak hydrophilic surfaces (silicon surfaces). This suggested that convective heat
transfer can be enhanced by increasing the surface hydrophilic capability.
A detailed 3-D analysis heat transfer in micro-channel heat sinks was provided by
Li et al. [39]. In their study the variation of the liquid thermophysical properties with
the bulk temperature was considered. For the 10mm long, 57?m wide, and 180?m deep
rectangular micro-channels analyzed, a variation in the reference temperature from 20?C
to 32?C changed the mean velocity from 1.11 to 1.31m/s. This resulted in a corresponding
change in Reynolds number from 96 to 144.
15
Ambatipudi and Rahman [2] investigated numerically the heat transfer in a ?101? silicon
substrate containing micro-channels. The effects of channel aspect ratio, Reynolds number,
and number of channels on the thermal performance of the device were investigated. It was
found that the Nusselt number is greater for a heat sink with a larger number of channels
and a larger Reynolds number. For Re = 673, the optimum channel depth that maximizes
Nusselt number occurred at 300?m.
Chein and Huang [5] studied the performance of micro-channel heat sinks using nanoflu-
ids as cooling fluids. The nanofluid was a mixture of pure water and nanoscale copper
particles with various volume fractions. It was found that nanofluids could enhance the
micro-channel heat sinks performance due to the increase in thermal conductivity of coolant
and the nanoparticle thermal dispersion effect.
Toh et al. [67] investigated the flow and heat transfer in micro-channels considering the
case of temperature dependent thermophysical properties. It was concluded that a drastic
change in the friction coefficient is due to the change in the viscosity of fluid. At lower
Reynolds numbers, the fluid attains higher temperatures leading to lower viscosities. This
results in lower pressure drop and hence the decrease in the value of the friction coefficient.
At higher Reynolds numbers, the fluid temperature at the exit is almost the same as the
inlet temperature. As a result the viscosity is almost constant and the friction coefficient
approaches the constant properties.
Two-phase micro-channel heat sinks (or micro heat pipes) are an alternative to single-
phase micro-channel heat sinks. They use the latent heat of vaporization to transfer heat
from the evaporator to the condenser, maintaining the sink at a uniform temperature.
Hassan et al. [18] present a state-of-the-art literature review of the research progress in the
16
field of micro-channel heat sinks. Earlier research works using single-phase coolants in their
heat sinks and more recent works using two-phase coolants were presented. They pointed
out that single-phase flow in micro-channel heat sinks requires high flow rates or smaller
hydraulic diameters, consequently resulting in larger pressure drops.
Groll et al. [15] provide an extended literature survey on the micro heat pipes used
for thermal control of electronic equipment. From a review of more than 100 papers micro
heat pipe sinks were classified and compared according to their performances, designs, and
thermophysical parameters. Common working fluids and wall materials were discussed and
possibilities for combinations were presented. As a measure of heat transport capacity a
liquid transport factor Nl was defined, Nl = ?hfg/?l, where ? is the surface tension, hfg
is the latent heat of vaporization, and ?l the kinematic viscosity of the liquid. This factor
indicated that for the temperature range of 0?C?100?C, typical for electronic applications
only acetone, methanol, ethanol, and water are suitable. Copper, aluminum, their alloys,
and stainless steel were mentioned as common heat pipewall materials. Tested combinations
of fluid/wall material with theoretically no degradation were referred to as compatible
combinations (e.g. water and either copper, stainless steel, nickel, titanium). Investigations
on the effect of fluid-solid contact angle were presented. It was concluded that the use of a
wetting fluid (small contact angle) should enhance the micro heat pipe performance. The
effect of the micro heat pipe dimensions was discussed in terms of Bond number. For the
micro heat pipe to operate properly regardless of orientation, the capillary forces should
overcome the gravitational forces leading to Bo < 1. Therefore micro heat pipes can not
be too long.
17
Experimental work on the performance of micro heat pipe arrays was conducted by
many researchers. Launay et al. [36] investigated the thermal behavior of two silicon micro
heat pipe arrays of triangular cross section. The first array made from two silicon wafers
included 55 triangular ethanol charged micro-channels, 230?m wide and 170?m deep. One
wafer contained the channels and the other one was used to seal them. The second array
made from three silicon wafers had 25 micro-channels 500?m wide and 342?m deep and
small triangular channels etched into the lower wafer which were used as arteries to drain the
liquid to the evaporator. Methanol was used as working fluid. With the assumption of no
heat losses and no heat conduction into the silicon wafer the effective thermal conductivity
of the second array was found to be 900W/mK, three times greater than the one of the
first array. This improvement was due to the increased flow cross sectional area and a
reduced pressure drop between the evaporator and condenser. The effect of liquid charge
on the performance of the first micro heat sink was also analyzed. It was found that
the enhancement factor for this micro heat pipe array, defined as the ratio between the
thermal conductivity of the charged and empty micro heat pipe, had a maximum value
of 7% corresponding to an optimum liquid charge of 24% of the channels? volume. This
enhancement factor was nearly constant for power inputs ranging from 0.5 to 5W. The
reason for this low improvement is that a large part of the heat is transferred by conduction
through the silicon.
Le Berre et al. [4] studied experimentally the performance of a micro heat pipe array
for various filling charges under various experimental conditions. The micro heat pipe
array was 20mm ? 20mm and consisted of 27 parallel triangular shaped channels, 500?m
wide and 350?m deep, giving a void fraction of 11%. They defined an effective thermal
18
conductivity as the thermal conductivity of a homogeneous material which would lead to
the same temperature field under the same conditions. The measurements indicated that
the use of micro heat pipes increases the effective thermal conductivity to about 200W/mK.
This represents an increase of about 67% of the thermal conductivity if compared to the
empty micro heat pipe array for which the thermal conductivity was determined to be
120W/mK. The results showed that the performance of the micro heat pipe array is
favored by decreasing the input heat flux or increasing the coolant temperature.
The research team from the Electrotechnics Laboratory of Grenoble performed exper-
imental and analytical work on flat rectangular micro heat pipes [13], [34], [27]. In their
design solution the evaporator (in contact with the electronic device) and the condenser (in
contact with a cold plate) were on opposite faces of the heat pipe. The proposed solution
was able to work in any position, even against gravity. They analyzed the behavior of the
heat spreader with and without working fluid. For the empty device the heat transfered
from the hot source to the cold plate was due only to conduction through the heat pipe
walls. The same experiment was run with the heat spreader filled with pure water and an
optimum liquid charge was found. The results showed that the thermal resistance between
the power device and the heat sink can be decreased by about 65% for the heat spreader
filled with 105?l of pure water compared to the thermal resistance obtained with the same
prototype without working fluid. The same team studied the advantages and disadvan-
tages of two-phase micro-channel heat sinks over single-phase micro-channel heat sinks [14].
They found that the pumping power required by the single-phase heat exchangers is 10
times greater than the power necessary for the operation of the two-phase counterparts.
The drawbacks of the two-phase heat exchangers were related to the working constraints
19
such as flow instabilities, local dry-out, critical heat flux and to the insertion of a condenser
in the cooling loop.
Suman and Hoda [63] analyzed the sensitivity of a V-shaped micro heat pipe to varia-
tions of the thermophysical properties and design parameters. They found that an increase
in the apex angle, viscosity, inclination, and length of the heat pipe reduces the performance
of the heat pipe. Suman et al [62] studied the performance of micro heat pipes having dif-
ferent polygonal cross sections. They showed that the critical heat input decreases with
increasing number of sides for regular polygons, due to an increase in the friction factor and
hence a decrease in liquid velocity and capillary pumping capacity. The radius of curvature
profile obtained by solving the coupled continuity, momentum, energy, and Laplace-Young
equations was used to predict the onset of the dry-out point and the propagation of the
dry-out length.
Microjet cooling devices are used as another solution for the thermal management
of microelectronic components. Kercher et al. [28] investigated the efficacy of synthetic
microjet technology for electronic cooling. Synthetic jets are formed from entrainment and
expulsion of the fluid in which they are embedded. The microjet cooling devices in this
study consist of a vibrating diaphragm, a diaphragm driver, and a housing with an orifice.
When the diaphragm is vibrated by the driver, air is drawn into the housing through the
orifice and ejected out from the same orifice, generating a series of vortex rings which
propagate away from the orifice. As a result a round turbulent jet is synthesized as these
vortex rings interact downstream. The cooling performance of the microjet cooling devices
was assessed using a thermal test die consisting of a diode-bridge temperature sensor and
a resistive heating element. It was shown that the performance of the microjet cooling
20
device was strongly dependent on the spacing between the orifice and the thermal die, the
maximum cooling performance being achieved for a spacing of 15mm ? 20mm. They also
concluded that the performance increases with increasing membrane resonance frequency
and increasing driving power. The influence of the orifice diameter on the performance was
also discussed and the 2.38mm orifice showed the best result, compared to the 1.98mm and
2.78mm orifices also used in this analysis. Fabbri and Dhir [11] performed a heat transfer
analysis using arrays of microjets which can improve the spatial uniformity of the heat
transfer coefficient. Ten different arrays were studied with the jet diameters varying from
69?m to 250?m and the Reynolds number from 73 to 3813. The best performance was
achieved using water jets of 173.6?m diameter and 3mm spacing, impinging at 12.5m/s on
a circular 19.3mm diameter copper surface.
A brief review of the results obtained by previous investigators on the effect of different
geometrical and thermophysical parameters of the heat sinks on their performance is pre-
sented in Tables 1.1-1.2. So far only experimental analysis has been done on the amount of
working fluid charging the MHP sinks. In this work an analytical methodology is devised
for assessing the performance of different configurations of MHP sinks and the amount of
liquid charge corresponding to a maximum heat transport capacity of the MHPs.
21
Author(s) Heat sink type/ Methodology/Results
Coolant
Hingorani et al. (1994) Planar spreader -optimum thickness (Bi and radius given)
-optimum radius (Bi and thickness given)
Peles et al. (2005) Micro pin fin heat sink -low thermal resistances
-flexibility in pin shape selection and spacing
Saha and Acharya (2003) Micro pin fin heat sink -increase in thermal performance when the
flow becomes unsteady
Tuckerman and Pease (1981) Micro-channel heat sink/ -optimum heat sink design (wfin/wch = 1,
water dissipating 790W/mm2, for ?T = 70?C and
?P = 30psi)
Knight et al. (1992) Micro-channel heat sink/ -allowed relaxation of the constraints used by
air Tuckerman and Pease obtaining designs with
improved performance
Weisberg et al. (1992) Micro-channel heat sink -optimum thermal resistance (wfin/wch = 1)
Vafai and Zhu (1999) 2-layer micro channel -reduced undesirable temperature gradients
heat sink in the streamwise direction
Ambatipudi and Rahman Micro-channel heat sink -Nu increases with increasing nr. of channels
(2000) -for Re = 673, maximum N is achieved for a
channel depth of 300?m
Ryu et al. (2002) Micro-channel heat sink -important thermal entrance effect when
the working fluid is water (Pr ? 7)
-thermal resistance increased by 15% when
nr. of channels was reduced by half
Toh et al. (2002) Micro-channel heat sink/ -investigated the temperature dependent
water thermophysical properties
Table 1.1: Literature review
22
Author(s) Heat sink type/ Methodology/Results
Coolant
Zhao and Lu (2002) Micro-channel heat sink -MHP modeled as lumped capacitance with
and effective thermal conductivity ? 10 times
that of Si or Cu)
Wu and Cheng (2003) Micro-channel heat sink/ -convective heat transfer enhanced by
water increasing the surface hydrophilic capability
Peng and Peterson (2004) Micro-channel heat sink -laminar flow: performance depends upon
aspect ratio H/W, wch/wfin
-turbulent flow: performance depends on
Z = min(H,W)/max(H,W), Zopt = 0.5
Chein and Huang (2005) Micro-channel heat sink/ -increase in thermal conductivity of coolant
nanofluids
Launay et al. (2004) MHP sink/ethanol, -1st array: 55 triangular channels,
methanol 230?m wide, 170?m deep
-2nd array: 25 triangular channels,
342?m deep, with arteries
-thermal conductivity of 2nd: 900W/mK
-optimum liquid charge for 1st: 24%
Suman and Hoda (2005) MHP sink/pentane -an increase in the apex angle, viscosity, tilt,
and length of MHP reduces the performance
Le Berre et al. (2006) MHP sink/methanol -tested a 20mm?20mm sink, with 27
triangular channels, 500?m wide, 350?m deep
-10% volume charge increases the effective
thermal conductivity by 67% over
the empty MHP array
Table 1.2: Literature review cont?d
23
Chapter 2
Optimal convective heat transfer coefficient of a planar spreader
2.1 Problem specification
Cooling of electronic devices requires the use of heat spreaders whose function is to
allow the spreading of the heat flux lines in the 3-D space and to increase the exchange area
with the coolant. This chapter is based on the work performed by Simionescu et al. [60]
in which the convective heat transfer coefficient on the bottom surface of a heat spreader
corresponding to a maximum heat removal from a heat source placed on the top surface of
the heat spreader is estimated (see Figure 2.1).
The elliptic partial differential equation for steady heat transfer is posed in a bounded,
3-D domain ? = ?c ??s of class C0,1 (with a Lipschitz continuous boundary; see [3]). The
temperature T satisfies
???(k?T) = qv in ? (2.1)
?T
?n
vextendsinglevextendsingle
vextendsinglevextendsingle
??\?
= 0 (2.2)
ks?T?n
vextendsinglevextendsingle
vextendsinglevextendsingle
?
= ? h(T|? ?T?) (2.3)
kc?Tc?n
vextendsinglevextendsingle
vextendsinglevextendsingle
S
= ? ks?Ts?n
vextendsinglevextendsingle
vextendsinglevextendsingle
S
(2.4)
Tc|S = Ts|S , (2.5)
24
Figure 2.1: Computational domain.
where k is the coefficient of thermal conductivity defined as:
k =
??
??
??
kc in ?c
ks in ?s.
Here ?c and ?s are the subdomains corresponding to the chip and spreader respectively, ??
is the boundary, n is the outward pointing normal vector to the boundary. The temperature
in the chip and spreader are Tc = T |?c and Ts = T |?s, respectively. For simplicity it is
assumed that the heat is generated uniformly by a source q within the chip and is removed
from the bottom surface of the heat sink (spreader) by convection, all the other surfaces
being insulated. Both the chip and the spreader are isotropic media having the thermal
conductivities kc and ks, respectively. The boundary conditions are of Robin type on ?,
the bottom surface of the spreader, corresponding to Newton?s law of cooling, Eq.(2.3),
where T? is the ambient temperature. The boundary conditions on the remaining bound-
ary, ??\? are of homogeneous Neumann type corresponding to an insulated boundary,
25
Eq.(2.2). Across the interface S between the two domains the temperature and heat flux
are continuous, Eq.(2.4)?(2.5).
Let H1 (?) be the usual Sobolev space of square integrable functions whose first dis-
tributional derivatives are also square integrable:
H1 (?) = {v:
integraldisplay
?
v2dx < ?,
integraldisplay
?
|?v|2 dx < ?}, (2.6)
equipped with the inner product
(u,v)H1(?) =
integraldisplay
?
uv +?u??vdx, (2.7)
and norm
bardblvbardblH1(?) = (v,v)1/2H1(?) =
parenleftbiggintegraldisplay
?
v2 +|?v|2 dx
parenrightbigg1/2
. (2.8)
We also make use of the space of square integrable functions L2(?):
L2(?) =
braceleftbigg
v ? ?:
integraldisplay
?
v2dx < ?
bracerightbigg
, (2.9)
equipped with the inner product
(u,v)L2(?) =
integraldisplay
?
uvdx, (2.10)
and norm
bardblvbardblL2(?) = (v,v)1/2L2(?) = (
integraldisplay
?
v2dx)1/2. (2.11)
26
Optimal control techniques have been applied successfully in the field of heat transfer
for solving boundarycontrol and parameter estimation problems. Cola?co and Orlande [7] es-
timated the boundaryheat fluxes at two surfaces of a cavity, by using simulated temperature
measurements taken from its interior. Martin and Dulikravich [43] estimated the convec-
tion coefficient using an inverse boundary element method. They used in their approach
internal temperature measurements or overspecified boundary conditions on the boundary
portion not exposed to convection. Silieti [58] used a genetic algorithm to estimate the
heat transfer coefficients. Li and Yan [38] analyzed the unsteady temperature distribution
in a turbulent forced convection between parallel plates taking the heat flux applied on
the upper wall as the control. Huang and Yeh [25] used the boundary heat flux to control
the entrance temperatures of concurrent flows between two adjacent parallel plate channels
of a concurrent flow heat exchanger. Another study by Huang and Yeh [24] presents the
temperature and moisture distributions inside a porous slab controlled by the convective
boundary conditions. The controls are the heat and mass transfer coefficients, respectively.
In the paper by Park and Lee [45] two different boundary control problems were analyzed
in which the solution of the heat equation was controlled by the boundary temperature and
the boundary heat flux, respectively. Lenhart and Wilson [37] studied the optimal control
of a heat transfer problem with convective boundary conditions. They focused their efforts
on establishing the existence and uniqueness of the solution of the optimality system. The
optimal control problem is solved using a conjugate gradient method which consists of the
following basic steps [23], [24], [25]:
1. Direct problem formulation
2. Optimal control problem formulation
27
3. Sensitivity problems formulation
4. Adjoint problems formulation
5. Gradient equations
6. Computational procedure
2.2 Existence and uniqueness of solution
We prove the existence and uniqueness of a weak solution to the problem (2.1)-(2.5)
by applying the Lax-Milgram lemma (see [3]).
Define H as:
H: = H1(?c ??s)
= {v: v|?c ? H1(?c), v|?s ? H1(?s) and ?(v|?c)
vextendsinglevextendsingle
vextendsingleS = ?(v|?s)
vextendsinglevextendsingle
vextendsingleS}
where ? is the trace operator (see [3]). To obtain a weak formulation, the differential
equation (2.1) is multiplied by v ? H and then integrated over ?. After applying Green?s
formula and the corresponding boundary conditions, we get:
integraldisplay
?
qvvdx =
integraldisplay
?
???(k?T)vdx (2.12)
= ?
integraldisplay
??
vk?T?nds+
integraldisplay
?
k?T ??vdx
= ?
integraldisplay
??\?
vk?T?nds?
integraldisplay
?
vks?T?nds+
integraldisplay
?
k?T ??vdx
=
integraldisplay
?
h(T ?T?)vds +
integraldisplay
?
k?T ??vdx.
28
Define the linear form:
L(v) =
integraldisplay
?
qvvdx +
integraldisplay
?
hT?vds (2.13)
and the bilinear form,
a(T,v) =
integraldisplay
?
hTvds +
integraldisplay
?
k?T ??vdx. (2.14)
For our particular problem qv = const for x ? ?c and qv = 0 for x ? ?s. Now the weak
form of problem (2.1)-(2.5) can be written as: Find T ? H such that
a(T,v) = L(v) ?v ? H. (2.15)
Obviously, L:H ? R is a continuous, linear form and a:H ?H ? R is a symmetric,
continuous, and coercive bilinear form, i.e. there exist positive constants K1, K2, K3 such
that:
|L(v)| ? K1bardblvbardblH(?) ?v ? H
|a(u,v)| ? K2 bardblubardblH(?) bardblvbardblH(?) ?u,v ? H
K3 bardblvbardbl2H(?) ? a(v,v) ?v ? H.
Hence the Lax-Milgram lemma guarantees the existence and uniqueness of a weak solution
to problem (2.1)-(2.5) (see [3]).
2.3 Optimal control for model problem
We now show the existence and uniqueness of an optimal control and describe a conju-
gate gradient algorithm which can be used to approximate solutions of the optimal control
29
problem. The conjugate gradient algorithm requires the solution of three coupled problems,
the direct problem, sensitivity problem and adjoint problem. We consider a model problem
with a rectangular chip and rectangular heat spreader (see Figure 2.1).
2.3.1 The direct problem
The governing heat equation and the corresponding boundary conditions (2.1)-(2.5) for
both the chip and the spreader, which are in thermal contact represent the direct problem.
The two domains are made from materials with different thermophysical properties. The
only heat source is within the chip.
The following dimensionless variables are introduced:
Tc ? Tc ?T?q
vD2s/ks
Ts ? Ts ?T?q
vD2s/ks
x ? xD
s
y ? yD
s
z ? zD
s
.
Dimensionless heat equations in the chip (?c) and the in the heat sink (?s) are respec-
tively:
?2Tc
?x2 +
?2Tc
?y2 +
?2Tc
?z2 = ?
ks
kc (2.16)
and
?2Ts
?x2 +
?2Ts
?y2 +
?2Ts
?z2 = 0 (2.17)
subject to the following boundary conditions:
?Tc
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls?lc2
= 0 ?Tc?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls+lc2
= 0 (2.18)
?Tc
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws?wc2
= 0 ?Tc?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws+wc2
= 0 (2.19)
30
?Ts
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ?Ts?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls
= 0 (2.20)
?Ts
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ?Ts?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws
= 0 (2.21)
?Tc
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1+dc
= 0 ?Ts?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=0
= BiTs|z=0 (2.22)
and the conditions of temperature and flux continuity at the interface:
Tc|z=1 = Ts|z=1 ?Ts?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
=
?
???
???
???
???
??
???
???
???
???
??
kc
ks
?Tc
?z
vextendsinglevextendsingle
vextendsinglez=1 if
?
???
???
???
???
Ls?lc
2 ? x ?
Ls+lc
2
Ws?wc
2 ? y ?
Ws+wc
2
0 otherwise.
(2.23)
Here Bi = hDs/ks is the Biot number which measures the ratio of conductive to convective
heat transfer resistance and represents the control variable.
The direct problem is solved using a finite difference scheme for an initial guess of the
control variable Bi.
2.3.2 The optimal control problem
In the optimal control problem, the optimal control Bi(x,y) is the unknown, and the
equations (2.16)-(2.23) represent constraints that have to be satisfied. The goal is to find
an optimal control that will maximize the amount of heat removed from the chip. This
31
statement is equivalent to minimizing the following functional:
J =
integraldisplay Ls+lc
2
L?l
2
integraldisplay Ws+wc
2
Ws?wc
2
integraldisplay 1+dc
1
T2c (x,y,z)dxdydz + ?2
integraldisplay Ls
0
integraldisplay Ws
0
Bi2(x,y)dxdy (2.24)
where Tc is the computed temperature within the chip subject to the constraints (2.16)-
(2.23) and ? is a weighting coefficient penalizing the control function. The goal of optimiza-
tion is to minimize the first term appearing in the definition of the functional J. The second
term is added to account for the cost of the control. The nonnegative penalty parameter ?
can be used to change the relative importance of the second term in Eq.(3.50). The squares
of the integrands guarantee the existence of a minimum. Note that ideally Tc should be 0,
in which case the chip temperature is equal to the ambient temperature.
2.3.3 Existence and uniqueness of an optimal control
In proving the existence and uniqueness of a solution to the optimal control problem
we follow the procedure presented by Lions [40], De los Reyes [42], and Gunzburger and
Lee [16].
Define Uad:= braceleftbigBi ? L2(?)bracerightbig as the set of admissible controls which is obviously not
empty. We want to prove the existence and uniqueness of a pair (T?,Bi?) ? H ?Uad that
minimizes the functional J, subject to equations (2.16)-(2.23).
Consider a minimizing sequence {(Tn,Bin)} such that
J(Tn,Bin) ?? inf
(T,Bi)?H?Uad
J(T,Bi) (2.25)
32
along with the weak formulation:
a(Tn,v) = L(v) ?v ? H(?). (2.26)
From the definition of J we have ?2 bardblBinbardbl2L2(?) ? J(Tn,Bin) < ?, which implies that Bin
is bounded. From the coercivity of the bilinear form a and continuity of the linear form L,
Tn is also bounded. So we may extract subsequences such that Bin ? Bi? weakly in L2(?),
Tn ? T? weakly in H(?), Tn ? T? strongly in L2(?) for some (T?,Bi?) ? H(?)?L2(?).
The last convergence follows from the compact embedding H(?) ???? L2(?). Now by the
weak lower semicontinuity of J we conclude that (T?,Bi?) is an optimal solution, i.e.
J(T?,Bi?) = inf
(T,Bi)?H?Uad
J(T,Bi). (2.27)
Thus we have shown that the optimal control exists. The uniqueness of the optimal
solution follows from the convexity of the functional J.
2.3.4 Conjugate gradient algorithm
A minimizing sequence is constructed using the conjugate gradient method and the
following iterative scheme (see [6]):
Bi(n+1) = Bi(n) ??(n)P(n) (2.28)
where ?(n) is the step size and P(n) is the search direction which is determined from:
P(n) = ?J(n) +?(n)P(n?1). (2.29)
33
The conjugate coefficient ?(n) is given by:
?(n) =
integraltextLs
0
integraltextWs
0 |?J
(n)(x,y)|2dxdy
integraltextLs
0
integraltextWs
0 |?J(n?1)(x,y)|2dxdy
with ?(0) = 0. (2.30)
The iterative process described by Eq. (2.28)?(2.30) requires the step size ?(n) and the
gradient of the functional, ?J(n). These parameters will be computed from the sensitivity
and adjoint problem, respectively.
2.3.5 Sensitivity problem and calculation of the search step size
When the control Bi(x,y) undergoes a small perturbation ?Bi(x,y), the temperatures
within the chip and spreader are varied by ?Tc(x,y,z) and ?Ts(x,y,z) respectively. If we
replace Tc in the direct problem (2.1)-(2.5) by Tc+?Tc, Ts by Ts+?Ts and Bi by Bi+?Bi
and then from the resulting equations subtract the equations describing the direct problem
the following sensitivity problem is obtained when the second order terms are neglected:
?2?Tc
?x2 +
?2?Tc
?y2 +
?2?Tc
?z2 = 0 (2.31)
?2?Ts
?x2 +
?2?Ts
?y2 +
?2?Ts
?z2 = 0 (2.32)
??Tc
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls?lc2
= 0 ??Tc?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls+lc2
= 0 (2.33)
??Tc
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws?wc2
= 0 ??Tc?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws+wc2
= 0 (2.34)
??Ts
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??Ts?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls
= 0 (2.35)
34
??Ts
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??Ts?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws
= 0 (2.36)
??Tc
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1+dc
= 0 ??Ts?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=0
= Bi?Ts|z=0 + ?BiTs|z=0 (2.37)
and the conditions of temperature and flux continuity at the interface:
?Tc|z=1 = ?Ts|z=1
??Ts
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
=
?
???
???
???
???
??
???
???
???
???
??
kc
ks
??Tc
?z
vextendsinglevextendsingle
vextendsinglez=1 if
?
???
???
???
???
Ls?lc
2 ? x ?
Ls+lc
2
Ws?wc
2 ? y ?
Ws+wc
2
0 otherwise
(2.38)
where ?Tc(x,y,z) and ?Ts(x,y,z) represent the sensitivity functions. The functional
J(n+1) for the (n+1)th iteration is written by replacing Bi(n+1) with the expression given
by Eq. (2.28):
J(n+1) =
integraldisplay Ls+lc
2
Ls?lc
2
integraldisplay Ws+wc
2
Ws?wc
2
integraldisplay 1+dc
1
T2c (x,y,z,Bi(n) ??(n)P(n))dxdydz
+ ?2
integraldisplay Ls
0
integraldisplay Ws
0
(Bi(n) ??(n)P(n))2dxdy. (2.39)
Here Bi(n+1) is a parameter. In the above and following equation we abuse notation by
adding a variable to Tc to emphasize its dependence on the control. Using a Taylor se-
ries expansion for Tc(x,y,z,Bi(n) ? ?(n)P(n)), the following expression is obtained for the
35
functional J at the (n+ 1)th iteration:
J(n+1) =
integraldisplay Ls+lc
2
Ls?lc
2
integraldisplay Ws+wc
2
Ws?wc
2
integraldisplay 1+dc
1
[Tc(x,y,z,Bi(n))??(n)?Tc(P(n))]2dxdydz
+ ?2
integraldisplay Ls
0
integraldisplay Ws
0
(Bi(n) ??(n)P(n))2dxdy. (2.40)
In the above equation Tc is the temperature distribution within the chip obtained by solving
the direct problem (2.16)-(2.23) using an initial guess of the control Bi, while ?Tc is one of
the sensitivity functions obtained from the solution of the sensitivity problem (2.31)-(2.38).
The search step size ?(n) is obtained from the line search at each iteration by taking the
derivative of Eq. (2.40) with respect to ?(n) and equating to 0:
?(n) =
integraltext Ls+lc2
Ls?lc
2
integraltext Ws+wc2
Ws?wc
2
integraltext1+ds
1 2Tc(Bi
(n))?Tc(P(n))dxdydz +?integraltextLs
0
integraltextWs
0 Bi
(n)P(n)dxdy
integraltext Ls+lc2
Ls?lc
2
integraltext Ws+wc2
Ws?wc
2
integraltext1+dc
1 2?T2c (P(n))dxdydz +?
integraltextLs
0
integraltextWs
0 (P(n))2dxdy
.
(2.41)
2.3.6 The adjoint problem
The adjoint problem is obtained by multiplying Eq. (2.16) and Eq. (2.17) by the
Lagrange multipliers ?1(x,y,z) and ?2(x,y,z), respectively (see [6]). The obtained equations
are integrated over the specified domains and the integrals are added to the functional J:
J =
integraldisplay Ls+lc
2
Ls?lc
2
integraldisplay Ws+wc
2
Ws?wc
2
integraldisplay 1+dc
1
T2c dxdydz + ?2
integraldisplay Ls
0
integraldisplay Ws
0
Bi2dxdy
+
integraldisplay Ls+lc
2
Ls?lc
2
integraldisplay Ws+wc
2
Ws?wc
2
integraldisplay 1+dc
1
?1(?
2Tc
?x2 +
?2Tc
?y2 +
?2Tc
?z2 +
ks
kc)dxdydz
+
integraldisplay Ls
0
integraldisplay Ws
0
integraldisplay 1
0
?2(?
2Ts
?x2 +
?2Ts
?y2 +
?2Ts
?z2 )dxdydz. (2.42)
36
To obtain the variation ?J we need to calculate J|Bi+?Bi ? J|Bi where J|Bi+?Bi is found
by replacing Tc by Tc + ?Tc, Ts by Ts + ?Ts and Bi by Bi + ?Bi in Eq. (2.42) and then
subtracting from the resulting expression the right hand side of Eq. (2.42):
?J|Bi =
integraldisplay Ls+lc
2
Ls?lc
2
integraldisplay Ws+wc
2
Ws?wc
2
integraldisplay 1+dc
1
2Tc?Tcdxdydz +?
integraldisplay Ls
0
integraldisplay Ws
0
Bi?Bidxdy
+
integraldisplay Ls+lc
2
Ls?lc
2
integraldisplay Ws+wc
2
Ws?wc
2
integraldisplay 1+dc
1
?1(?
2?Tc
?x2 +
?2?Tc
?y2 +
?2?Tc
?z2 +
ks
kc)dxdydz
+
integraldisplay Ls
0
integraldisplay Ws
0
integraldisplay 1
0
?2(?
2?Ts
?x2 +
?2?Ts
?y2 +
?2?Ts
?z2 )dxdydz. (2.43)
Integrating Eq.(2.43) by parts and applying the boundary conditions of the sensitivity
problem the following adjoint problem is obtained:
?2?1
?x2 +
?2?1
?y2 +
?2?1
?z2 + 2Tc = 0 (2.44)
?2?2
?x2 +
?2?2
?y2 +
?2?2
?z2 = 0 (2.45)
??1
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls?lc2
= 0 ??1?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=Ls+lc2
= 0 (2.46)
??1
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws?wc2
= 0 ??1?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=Ws+wc2
= 0 (2.47)
??2
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??2?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L
= 0 (2.48)
??2
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??2?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=W
= 0 (2.49)
??1
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1+dc
= 0 ??2?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=0
= Bi?2|z=0 (2.50)
37
?1|z=1 = kck
s
?2|z=1 ??2?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
=
?
???
???
???
???
??
???
???
???
???
??
??1
?z
vextendsinglevextendsingle
vextendsinglez=1 if
?
???
???
???
???
Ls?lc
2 ? x ?
Ls+lc
2
Ws?wc
2 ? y ?
Ws+wc
2
0 otherwise.
(2.51)
The following integral is left after performing the integration by parts in Eq. (2.43)
and identifying the adjoint problem:
?J|Bi =
integraldisplay Ls
0
integraldisplay Ws
0
(?Bi??2(x,y,0)Ts(x,y,0))?Bidxdy. (2.52)
By the definition given in [1], the functional increment at Bi(x,y) can be written as:
?J|Bi =
integraldisplay Ls
0
integraldisplay Ws
0
?J?Bidxdy. (2.53)
Hence we obtain the following expression for the gradient:
?J = ?Bi??2(x,y,0)Ts(x,y,0). (2.54)
2.3.7 Computational algorithm
The computational procedure for the solution of this problem may be summarized as
follows:
step 1 Solve the direct problem given by Eq. (2.16)?(2.23) to obtain Tc and Ts for an initial
guess of the control Bi;
38
step 2 Solve the adjoint problem given by Eq. (2.44)?(2.51) to obtain ?1 and ?2;
step 3 Compute the gradient of the functional ?J from Eq. (2.54);
step 4 Compute the conjugate coefficient ?(n) from Eq. (2.30) and the search direction P(n)
from Eq. (2.29);
step 5 Set ?Bi(x,y) = P(n)(x,y) and solve the sensitivity problem given by Eq. (2.31)?
(2.38) to obtain ?Tc and ?Ts respectively;
step 6 Compute the search step size ?(n) from Eq. (2.41);
step 7 Compute the new estimate for Bi from Eq. (2.28) and go back to step 1 until a
stopping criterion is achieved. The stopping criterion for our numerical experiments
was chosen as
vextenddoublevextenddouble
vextenddoubleBi(n+1) ?Bi(n)
vextenddoublevextenddouble
vextenddouble? ? , where ? is a prescribed tolerance which was set
in our experiments to 10?6.
2.4 Implementation and numerical procedure
A finite difference scheme was implemented to solve the direct, sensitivity, and adjoint
problems. The centered approximation for the derivatives was used leading to the following
discrete heat transfer equations for the direct problem:
Ti?1,j,kc ?2Ti,j,kc +Ti+1,j,kc
?x2 +
Ti,j?1,kc ?2Ti,j,kc +Ti,j+1,kc
?y2
+ T
i,j,k?1c ?2Ti,j,kc +Ti,j,k+1c
?z2 = ?
ks
kc
39
Ti?1,j,ks ?2Ti,j,ks +Ti+1,j,ks
?x2 +
Ti,j?1,ks ?2Ti,j,ks +Ti,j+1,ks
?y2
+ T
i,j,k?1s ?2Ti,j,ks +Ti,j,k+1s
?z2 = 0
at x = 0 T0,j,ks = T2,j,ks
at x = Ls Tm?1,j,ks = Tm+1,j,ks
at y = 0 Ti,0,ks = Ti,2,ks
at y = Ws Ti,n?1,ks = Ti,n+1,ks
at z = 0 Ti,j,0s ?Ti,j,2s2?z = BiTi,j,1s
at z = 1 Ti,j,ps = Ti,j,pc
Ti,j,p?1s ?Ti,j,p+1s = kcks(Ti,j,p?1c ?Ti,j,p+1c )
at x = Ls?lc2 T
m?m0
2 ,j,kc = T
m?m0
2 +2,j,kc
at x = Ls+lc2 T
m+m0
2 ?1,j,kc = T
m+m0
2 +1,j,kc
at x = Ws?wc2 Ti,
n?n0
2 ,kc = Ti,
n?n0
2 +2,kc
at x = Ws+wc2 Ti,
n+n0
2 ?1,kc = Ti,
n+n0
2 +1,kc
at z = 1 +dc Ti,j,p+p0?1c = Ti,j,p+p0+1c
where Ti,j,ks and Ti,j,kc represent the spreader and chip temperature, respectively correspond-
ing to the (i,j,k) location of the mesh, m, n, p, represent the number of grid points in x,
y, and z direction in the spreader and m0, n0, p0, represent the number of grid points in
x, y, and z direction in the chip. The step sizes for the two domains have the same values.
40
0 1000 2000 3000 4000 5000
0
500
1000
1500
2000
2500
3000
3500
4000
4500
5000
nz = 32699
Figure 2.2: Structure of matrix A.
The discretization yields a linear system of equations:
AT = b.
Here A is a d ? d matrix having the banded structure presented in Figure 2.2, where
d = (m + 1)(n + 1)(p + 1) + (m0 + 1)(n0 + 1)p0, T is the d vector of the temperatures at
the grid points and b is the d vector of the free terms corresponding to the heat source
component. Similar discretizations were used for the sensitivity and adjoint problems. A
Matlab code was implemented to solve the three coupled problems, using a mesh of 31?31?5
points for the spreader and 11?11?4 points for the chip.
2.5 Results and discussions
We now present some simulation results. Parameter values that were used in the
simulation were:
? volumetric heat source: qv = 144 W/mm3
41
? heat source (chip) dimensions: 3mm?3mm?0.15mm
? heat sink (spreader) dimensions: 9mm?9mm?0.3mm
? chip thermal conductivity: kc = 0.1 W/mmK
? spreader thermal conductivity: ks = 0.1412 W/mmK
These values are typical of the ones for applications we have in mind. The controlled and
uncontrolled solutions for the temperature distributions are compared for different values of
?. The uncontrolled solution is obtained by solving the direct problem presented in section
2.3.1 with a constant Biot number over the entire surface ?. The value used in each case
(for each value of the parameter ?) would give the same cost as the cost of the optimal
Bi. The distribution of the temperature difference T ? T? at the interface between the
chip and spreader for the controlled case is lower and displays a flatter profile than the
uncontrolled solution as we can see in Figure 4.1. We see that the difference between the
two distributions increases with increasing penalty parameter ? (increasing the cost of the
control). Small variations of the penalty parameter ? introduce significant variations of the
temperature distribution inside the domain. The greater the value assigned to ? the larger
the cost of employing a control, resulting in a smaller Bi hence a smaller contribution of
convective cooling.
Figure 4.2 shows distributions of the Biot number for different values of the penalty pa-
rameter, ? = 0.1,0.5,1,5,50. The corresponding distributions of the temperature difference
at the chip-spreader interface are presented in Figure 4.3. The profile of Bi resembles the
temperature profile at any horizontal cross section through the spreader. The smaller the
penalty parameter ?, the higher the Biot number and the lower the domain temperature.
42
The temperature of the hottest point, which is the center of the chip?s upper face,
corresponding to ? = 50 is approximately 50K higher than the one corresponding to ? =
0.1. These temperatures are approximately 150K and 50K lower than the uncontrolled
temperatures, respectively.
Contour plots of the temperature distributions, trough a vertical symmetry plane are
presented in Figure 4.4 and show that the temperature difference between the hottest and
coolest point increases with increasing ? and the contour curves are flatter for small values
of ?.
43
Chapter 3
Optimal design of a micro heat pipe
3.1 Problem specification
In this chapter the analysis is extended to heat sinks equipped with micro heat pipes.
The objective of this analysis is to estimate the geometry of a micro heat pipe (MHP) and
the amount of cooling fluid inside the pipes that will maximize the heat removed from a heat
source placed on the top surface of the heat sink (see Figure 3.1). Micro heat pipes are small
sealed devices filled with a working fluid whose phase change is used to transfer thermal
energy. A single micro heat pipe consists of a small non circular channel that uses the
sharp angled corners as liquid arteries. The vaporization-condensation process causes the
liquid-vapor interface to change along the pipe, resulting in a capillary pressure difference.
The capillary pumping pressure drives the liquid from the condenser back to evaporator.
An optimal control technique similar to that presented in chapter 2 is used to solve
two problems. In the first problem the solution of the heat equation is controlled by the
convective boundary condition taking the Biot number and the geometrical parameters
of the channel as controls. In the second problem an optimal Bi number distribution
corresponding to the MHP maximum heat transport capacity is found. The optimal heat
transfer coefficient can be obtained by prescribing the height and width of the channels,
the number of channels, and the amount of working fluid that charges the channels. The
importance of the channel geometry on the performance of a heat sink was pointed out by
many previous investigators as presented in chapter 1. Only experimental work was done
to study the influence of liquid charge on the performance of MHP sinks [36] and it was
44
Figure 3.1: Schematic of a MHP sink.
shown that a too large amount of fluid leads to a condenser flooding, while a too low fluid
charge leads to an evaporator dry-out and an increase in channel wall temperature.
3.2 The direct problem
3.2.1 Heat transfer equations corresponding to the GaAs slab and silicon wafer
In this analysis the silicon wafer dimensions are constant and all the channels have
a uniform rectangular cross section of width 2l and height H1. A constant heat flux q is
supplied to a GaAs substrate through a Lt ?w2 heating area. The GaAs slab of width Lt,
length w1, and height H3 is placed on the top of the silicon wafer, covering the evaporator
area of the micro heat pipes. Taking advantage of symmetry we consider the analysis of a
cell consisting of half of the channel and half of the separating fin, see Figure 3.2.
45
Figure 3.2: Computational domain.
The heat transfer equations and the corresponding boundary conditions are written
for each subdomain depicted in Figure 3.2, i.e. GaAs slab, spreader, and fin. It is assumed
that the heat entering through the Lt?w2 heating area it is removed only from the channel
surface through the evaporator area of the MHP, all the other surfaces being insulated. The
dimensionless mathematical formulation of the heat conduction through the GaAs slab,
silicon substrate, and fin is obtained by introducing the following dimensionless variables:
TA ? TA ?TfqH
2/kSi
TB ? TB ?TfqH
2/kSi
TC ? TC ?TfqH
2/kSi
x ? xH
2
y ? yH
2
z ? zH
2
.
Here TA,TB,TC,Tf are the temperatures corresponding to the GaAs slab, silicon substrate,
fin, and cooling fluid respectively, q is the heat flux entering the heat sink, kSi is the
coefficient of thermal conductivity for silicon, and the reference length H2 is the thickness
of the silicon wafer. In what follows it is understood that all quantities - unless specified
46
otherwise - are in dimensionless form. The governing equations are:
?2TA
?x2 +
?2TA
?y2 +
?2TA
?z2 = 0 (3.1)
?2TB
?x2 +
?2TB
?y2 +
?2TB
?z2 = 0 (3.2)
?2TC
?x2 +
?2TC
?y2 +
?2TC
?z2 = 0 (3.3)
subject to the following boundary conditions:
?TA
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ?TA?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L
= 0 (3.4)
?TA
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ?TA?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=w1
= 0 (3.5)
?TB
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ?TB?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?l
= 0 (3.6)
?TB
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ?TB?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=W
= 0 (3.7)
?TC
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ?TC?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?l
=
?
???
???
?Bi TC|x=L?l if 0 ? y ? Le
0 if Le < y ? W
(3.8)
?TC
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ?TC?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=W
= 0 (3.9)
?TC
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=0
= 0 ?TA?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1+H3
=
?
???
???
kSi
kGaAs if 0 ? y ? w2
0 if w2 < y ? w1
(3.10)
47
and the conditions of temperature and flux continuity across the GaAs slab-silicon substrate
interface,
TB|z=1 = TA|z=1 if 0 ? y ? w1 (3.11)
?TB
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
=
??
??
???
kGaAs
kSi
?TA
?z
vextendsinglevextendsingle
vextendsinglez=1 if 0 ? y ? w1
0 if w1 < y ? W
(3.12)
and across the silicon substrate-fin interface,
TB|z=H1 = TC|z=H1 if 0 ? x ? L?l (3.13)
?TB
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=H1
=
?
???
???
???
???
??
???
???
???
???
??
?TC
?z
vextendsinglevextendsingle
vextendsinglez=H
1
if 0 ? x ? L?l
BiTB|z=H1 if
?
??
???
L?l < x ? L
0 ? y ? Le
0 if
?
???
???
L?l < x ? L
Le < y ? W.
(3.14)
For the convective boundary conditions corresponding to the silicon substrate and fin it is
assumed that the convective heat transfer coefficients are equal, h1 = h2 = h, and vary
only in y-direction. This is a simplifying assumption. In reality h is not constant at any
cross section of the channel, it depends on the channel wall temperature which has a 2-
D distribution. The convective boundary conditions are expressed in terms of the Biot
number:
Bi = hH2k
Si
. (3.15)
48
IftheBiot numberisconstant along the micro-channel crosssection, ? = {x : L?l ? x ? L}
uniontext{z : 0 ? z ? H
1}, and Tw is the channel wall temperature the heat removed from the
source is:
integraldisplay W
0
integraldisplay
?
qwdsdy =
integraldisplay W
0
integraldisplay
?
?Tw
?n dsdy
=
integraldisplay W
0
parenleftBiggintegraldisplay H
1
0
? ?TC?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?l
dz +
integraldisplay L
L?l
?TB
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
x=H1
dx
parenrightBigg
dy (3.16)
=
integraldisplay W
0
Bi(y)
parenleftBiggintegraldisplay H
1
0
TC(y,z)|x=L?l dz +
integraldisplay L
L?l
TB(x,y)|z=H1 dx
parenrightBigg
dy.
At any cross section y, the heat flux is:
qw|y = 1?
integraldisplay
?
qwds = Bi(y)l +H
1
(
integraldisplay H1
0
TC(y,z)|x=L?l dz +
integraldisplay L
L?l
TB(x,y)|z=H1 dx). (3.17)
Hence the local Biot number is:
Bi(y) = qw|y (l +H1)integraltextH1
0 TC(y,z)|x=L?l dz +
integraltextL
L?l TB(x,y)|z=H1 dx
. (3.18)
3.2.2 1-D Micro heat pipe model
The micro heat pipe, which was first proposed by Cotter [8], has the advantage of large
latent energies associated with phase change. Another advantage of the micro heat pipes is
a nearly uniform temperature maintained throughout the device due to the phase change.
An overview of the state-of-the-art micro heat pipe sinks fabrication and analysis is provided
by Tien et al. [66] and Faghri [12]. Numerous analyses have been performed on micro heat
pipes, considering the thermodynamics of the liquid-vapor interface [61], [64], [17], [31], the
Marangoni effects [64], and the disjoining pressure[29], [64], [57]. Sartre et al. [57] developed
49
a 3-D steady state model for predicting the heat transfer in a micro heat pipe array. They
assumed that the extended meniscus is typically divided into three regions: a macroscopic
region with a meniscus of constant mean curvature which is dominated by the capillary
forces (intrinsic meniscus region), a transition region or microregion (evaporating thin film
region) and a microscopic region which consists of an adsorbed film (non-evaporating region)
where the disjoining forces predominate. The results presented for an aluminum/ammonia
triangular micro heat pipe array show that the major part of the total heat input goes
through the microregion. The mathematical model presented by Khrustalev and Faghri
[29], [30] accounted for the effects of interfacial thermal resistance, disjoining pressure,
and surface roughness for a given meniscus contact angle. The free surface temperature
of the liquid film was determined using the extended Kelvin equation and the expression
for interfacial resistance given by the kinetic theory. In this study a 1-D model similar
to that proposed by Longtin et al. [41] for triangular cross section micro heat pipes is
employed. Longtin?s model used a uniforminput heat flux while in this analysis the heat flux
has a 1-D distribution. The model is developed for the evaporator and adiabatic sections
and solved to yield pressure, velocity, and film thickness information in the lengthwise
direction of the pipe. Equations for the conservation of mass, momentum, and energy are
developed considering only axial variations. During operation the working fluid recedes into
the corners of the pipe, generating the necessary capillary driving pressure. Evaporation
in the evaporator and condensation in the condenser cause the liquid to accumulate in the
condenser section of the device (Figure 3.3).
The amount of liquid- hence the interfacial radius of curvature- varies as a function
of axial position y. The interfacial radius of curvature is related to the pressure difference
50
Figure 3.3: Schematic of a MHP.
between the liquid and the vapor phases by the Laplace-Young equation,
Pv(y)?Pl(y) = ?r(y) (3.19)
or by differentiating with respect to y,
dPv
dy ?
dPl
dy = ?
?
r2
dr
dy. (3.20)
Thederivation of theconservation of mass and momentum equations is performedemploying
the following assumptions:
? Both liquid and vapor flows are incompressible and the Reynolds number for both
flows does not exceed 50;
? Steady state operation;
? Constant fluid properties;
51
? Negligible viscous dissipation as a consequence of the small velocities in both the
liquid and the vapor phases;
? Constant vapor temperature Tv;
? The mean radius of curvature is approximately equal to the radius of curvature normal
to the pipe axis;
? The interfacial radius of curvature is constant at any given y location.
Figure 3.4: Control volume used to derive the conservation equations
The heat pipe is divided into a series of small control volumes of length dy for which
the conservation principles are applied. Some geometrical parameters necessary in the
derivation of the governing equations are defined as functions of the contact angle ?, pipe
corner angle ? = ?/2 and meniscus radius r(y) (see Figure 3.5).
The liquid cross sectional area (for the 2 corners of our computational domain) is:
Al = 2(A? ?Aa) = 2(?
2
2 ?Aa)
= 2(?
2
2 ?Asect +
r2sin?
2 )
= r2(2sin2?2 +sin? ??) = clr2(y). (3.21)
52
Figure 3.5: A cross section of the MHP with liquid recessed into the corner
The vapor cross sectional area is:
Av = lH1 ?Al = lH1 ?clr2(y). (3.22)
The wall area in contact with the liquid phase and vapor phase, respectively for a control
volume of thickness dy is given by:
Alw = 4?dy = 4?2sin?2rdy = clwrdy (3.23)
Avw = (2l +H1 ?4?)dy = (2l +H1 ?4?2sin?2r)dy = (2l +H1 ?clwr)dy (3.24)
and the area of the liquid-vapor interface for a control volume dy is given by:
Ai = 2r?dy = cirdy. (3.25)
53
In the above equations cl, clw, and ci are geometrical constants defined by:
cl = (2sin2?2 +sin? ??), clw = 4?2sin?2, ci = 2?. (3.26)
It should be mentioned that the radius of curvature must increase monotonically from the
evaporator to the condenser to avoid sign changes in the pressure gradient. A continuously
increasing radius of curvature results in a continuously increasing liquid cross sectional area.
If two liquid films meet in the evaporator or adiabatic section, the device will fail to operate
because dr/dy changes sign. The radius of curvature at this point is called the maximum
radius of curvature [41] and is given by:
rmax = min(l,H1/2)?2sinparenleftbigpi
4 ??
parenrightbig. (3.27)
Conservation of energy. The heat entering the computational domain through the
L ? w2 area is removed by the heat pipe through the evaporator section. It is assumed
that the heat load is uniformly distributed between the four corners of the pipe and is
entirely removed through the evaporator section. The following equation is obtained when
applying an energy balance on the computational domain (see Figure 3.2):
qLw2 =
integraldisplay Le
0
integraldisplay
?
qwdsdy
=
integraldisplay Le
0
h(y)
integraldisplay
?
(Tw ?Tf)dsdy (3.28)
where q is the heat flux applied to the heat sink, qw is the heat flux at channel wall, Tw
is the channel wall temperature, h is the convective heat transfer coefficient, and Le is
54
the evaporator length. Since the micro heat pipe is nearly isothermal during operation,
the liquid film is very thin, and the Reynolds numbers are very low we can neglect the
conduction, convection and viscous dissipation in the liquid. The heat is transferred only by
vaporization at the liquid-vapor interface. The assumption of constant vapor temperature
along the pipes implies that only the liquid phase has an expression for the conservation of
energy. Applying an energy balance on a control volume of length dy, we get:
qw|y dy = h(y)
parenleftBiggintegraldisplay L
L?l
parenleftBig
TB|z=H1 ?Tf
parenrightBig
dx+
integraldisplay H1
0
parenleftBig
TC|x=L?l ?Tf
parenrightBig
dz
parenrightBigg
dy
= qH22Bi(y)
parenleftBiggintegraldisplay H
1
0
TC
vextendsinglevextendsingle
vextendsinglex=L?l dz +
integraldisplay L
L?l
TB
vextendsinglevextendsingle
vextendsinglez=H
1
dx
parenrightBigg
dy
= ?lVilAihfg = ?lVilhfgcirdy (3.29)
where ?l is the liquid density, Vil is the liquid interface velocity, and hfg is the latent heat
of vaporization. The overscored quantities are in dimensionless form. From the equation
of energy conservation we can obtain the expression for the interfacial liquid velocity (in
dimensional form):
Vil =
qBi(y)
parenleftbiggintegraltext
H1
0 TC(?y,?z)
vextendsinglevextendsingle
vextendsinglex=L?l dz +integraltextLL?l TB(?x, ?y)
vextendsinglevextendsingle
vextendsinglez=H
1
dx
parenrightbigg
?lhfgci?r (3.30)
Conservation of mass. Referring to Figure 3.6 and using the geometrical parameters
derived above, Eq.(3.21)-(3.26), after performing a mass balance on a control volume of
thickness dy the following continuity equations are obtained for the liquid and vapor, re-
spectively:
55
Figure 3.6: Conservation of mass
2Uldrdy +rdUldy ? cic
l
Vil = 0 (3.31)
parenleftBig
lH1 ?clr2
parenrightBigdUv
dy ?2clr
dr
dyUv ?
?l
?vcirVil = 0 (3.32)
The continuity of mass across the interface, ?lVil = ?vViv is used in Eq.(3.32) to relate
the vapor interface velocity to that of the liquid.
Conservation of momentum. Due to the small velocities in both phases the convective
terms in the momentum equations are neglected. The body force (gravity) is significant
only for the liquid phase. The surface forces on the control volume are the pressures acting
on the cross sectional area and the shear stresses encountered at the liquid-vapor interface
and wall (see Figure 3.7). A balance of the body and surface forces on the control volume
56
Figure 3.7: Conservation of momentum
yields the following equations for the liquid and vapor respectively:
2Pldrdy + dPldy r ? cic
l
?li ? clwc
l
?lw +?lgrsin?1 = 0 (3.33)
?2Pvclrdrdy +
parenleftBig
lH1 ?clr2
parenrightBigdPv
dy + (2l +H1 ?clr)?vw +cir?vi = 0. (3.34)
In the above equations ?lw and ?vw are the liquid-wall and vapor-wall shear stresses, ?vi
is the shear stress at the liquid-vapor interface, and ?1 is the pipe inclination angle with
respect to gravity. Expressions for the shear stresses are derived by assuming that both
flows are fully developed. This assumption is justified since the convective terms are small
and the liquid and vapor cross sections vary slowly with y. The ratio of the vapor velocity
to the liquid velocity scales with the inverse ratio of the densities, ?l/?v. For this reason,
57
from the perspective of the vapor the liquid can be treated as another section of the wall.
This assumption was also employed by Thomas et al. [65] in their study of the two phase
laminar flow in trapezoidal grooves and by Khrustalev and Faghri [31] who modeled the
opposite flows of the vapor and liquid in miniature heat pipes. The interfacial shear stress
is computed for the vapor by assuming the liquid to be stationary. Since the liquid and
vapor shear stresses at the interface are equal and have opposite directions, the liquid shear
follows immediately. The interfacial shear stress due to countercurrent flow in a heat pipe
decreases the maximum heat transport. For cases in which the vapor velocity is high this
effect is more pronounced. Therefore the amount of liquid charge has a significant impact
on the heat pipe performance. The hydraulic diameters for the two phases are computed
using the parameters derived earlier in Eq.(3.21)-(3.26):
DHl = 4Al4? = rcl?2sinparenleftbigpi
4 ??
parenrightbig (3.35)
DHv = 4Av2l +H
1 ?clwr +Ai/dy
= 4
parenleftbiglH
1 ?clr2
parenrightbig
2l +H1 ?clwr +cir. (3.36)
The shear stresses are expressed as:
?vi = ?vw = 12?vU2vfv, fv = kvRe
v
, Rev = ?vUvDHv?
v
(3.37)
?li = ??vi (3.38)
?lw = 12?lU2l fl, fl = klRe
l
, Rel = ?lUlDHl?
l
(3.39)
58
The component of the momentum due to the mass entering or leaving the interface during
the phase change is neglected [47] . The constants kv and kl depend upon the geometry of
the ducts. For the liquid flow a rectangular geometry is assumed. The following relation
was obtained by Wu and Cheng [73] for a trapezoidal cross section:
k = 11.43 + 0.8e2.67Wb/Wt (3.40)
where Wb and Wt are the bottom and top widths of the trapezoidal profile, with Wb ? Wt.
If in the above equation we make Wb = Wt, corresponding to a rectangular profile a value
of approximately 13 is obtained for kl. For the vapor flow the constant kv is obtained by
averaging between the rectangular cross section displayed in the evaporator section and the
almost circular profile attained at the end of the adiabatic section. For the circular cross
section k = 16 and hence the value resulting for kv is 14.5. Replacing the expressions for ?lw,
?li, ?vw, and ?vi into the momentum Eq.(3.33)-(3.34) the following equations are obtained:
2Pldrdy + dPldy r ? cic
l
kvUv?v
DHv ?
clw
cl
klUl?l
DHl +?lgrsin?1 = 0 (3.41)
?2Pvclrdrdy +
parenleftBig
lH1 ?clr2
parenrightBigdPv
dy + (2l +H1 ?clwr +cir)
kvUv?v
DHv = 0 (3.42)
The transport Eq.(3.20), (3.31), (3.32), (3.33), (3.34) corresponding to the micro heat
pipe are cast into dimensionless form by introducing the following dimensionless variables:
?y ? yH
2
, ?r ? rH
2
, ?DH ? DHH
2
?Ul ? Ul
Uref =
Ul
Q/parenleftbig?lD2Hhfgparenrightbig,
?Uv ? Uv
Uref =
Uv
Q/parenleftbig?lD2Hhfgparenrightbig,
59
?Vil ? Vil
Uref =
Bi(?y)
parenleftBigintegraltext ?H
1
0 ?TC(?y, ?z)
vextendsinglevextendsingle
?x=?L??l d?z +
integraltext ?L
?L??l ?TB(?x, ?y)
vextendsinglevextendsingle
?z= ?H1 d?x
parenrightBig
ci?r
D2H
Lw2 (3.43)
?Pl ? Pl
?/DH =
Pl
Pref
?Pv ? Pv
?/DH =
Pv
Pref .
whereDH = 2lH1/(l +H1) is the hydraulic diameter of the empty pipe, Le isthe evaporator
length, Q is the input power, Uref = Q/parenleftbig?lD2Hhfgparenrightbig is the reference velocity, and Pref =
?/DH is the reference pressure. The equations in dimensionless form for the micro heat
pipe are (the overbars are omitted):
dPv
dy ?
dPl
dy =
2lH1
l +H1
1
r2
dr
dy (3.44)
2Uldrdy +rdUldy ? cic
l
Vil = 0 (3.45)
parenleftBig
lH1 ?clr2
parenrightBigdUv
dy ?2clr
dr
dyUv ?
?l
?vcirVil = 0 (3.46)
2Pldrdy + dPldy r ? Cacic
l
kv?v?
l
2l +H1 ?clwr +cir
4(lH1 ?clr2) UvDH
? Caclwc2
l
kl
?2sin?
2
r UlDH +Bo
DHrsin?
W2 = 0 (3.47)
?2Pvclrdrdy +
parenleftBig
lH1 ?clr2
parenrightBigdPv
dy +Ca
(2l +H1 ?clwr +cir)2
4(lH1 ?clr2) kvDHUv = 0. (3.48)
In the above equations (3.44)-(3.48), the geometrical parameters of the channel l, H1, W,
DH are in dimensionless form. The dimensionless groups introduced,
Ca = ?lUref? , Bo = ?lgW
2
?
60
represent the capillary number and Bond number, respectively. The capillary number Ca
is the ratio of the viscous force to surface tension force and the Bond number Bo measures
the contribution of the gravitational and capillary forces. Since the transport equations for
the heat pipe are first order ordinary differential equations, boundary conditions are needed
at only one point. The solution of these equations has two parts, one corresponding to
the evaporator section and one to the adiabatic section. The boundary conditions for the
adiabatic section are taken from the solution of the evaporator section at the evaporator-
adiabatic interface. The boundary conditions used at the end of the evaporator (y = 0) in
nondimensional form are:
r|y=0 = r0, Ul|y=0 = 0, Uv|y=0 = 0, (3.49)
Pv|y=0 = Pv(Tv), Pl|y=0 = Pv|y=0 ? DHr
0
.
The boundary vapor pressure Pv|y=0 is taken to be the saturation pressure of the vapor at
temperature Tv. The liquid boundary pressure is related to the vapor boundary pressure
by the Laplace-Young equation. The value for the meniscus radius at y = 0, r0 is given as
an initial guess [41].
The direct problem associated with the mathematical formulation given by Eq.(3.1)-
(3.14), involves the determination of temperature fields in the GaAs substrate, silicon sub-
strate, and fin from the knowledge of domain geometry and the physical properties of GaAs
and silicon.
61
3.3 The optimal control problem with Bi and geometrical parameters as con-
trols
An optimum Bi distribution for a finned heat sink was found previously by Simionescu
and Harris [59]. The analysis presented here uses two controls, i.e. the Biot number and
the aspect ratio, being an extension of the previous work. In the optimal control problem
Bi and p = l/H1 are the unknowns, where p is the channel aspect ratio, and Eq.(3.1)-(3.14)
are constraints that have to be satisfied. The problem is to find the optimum channel
geometry and the optimum Biot number that maximize the amount of heat removed from
the heat source (or minimize the thermal resistance). The difference between this problem
and the classic fin shape optimization problem studied previously [74], [22], [52] consists
in the consideration of the heat removed from the interfin area. This amount of heat was
neglected in previous works, but it is taken into account in this study since we made the
assumption of uniform heat load distribution between the four corners of the pipe. The
following functional is considered for minimization in our optimal control problem:
J =
integraldisplay W
0
integraldisplay L
0
integraldisplay 1+H3
1
T2Adxdydz + ?2
integraldisplay Le
0
Bi2dy. (3.50)
Here TA is the computed temperature within the GaAs slab, subject to the constraints
dictated by the direct problem (3.1)-(3.14), and ? is a weighting coefficient penalizing Bi.
The goal of optimization is to minimize the first term appearing in the definition of the
functional J. The second term is added to account for the cost of the control parameter
Bi. The nonnegative penalty parameter ? can be used to change the relative importance
of the last term in Eq.(3.50).
62
3.3.1 Conjugate gradient algorithm
Minimizing sequences are constructed using the conjugate gradient method and the
following iterative scheme (see [6]):
p(n+1) = p(n) ??(n)1 P(n)1 (3.51)
Bi(n+1) = Bi(n) ??(n)2 P(n)2 (3.52)
where ?(n)1 and ?(n)2 are the step sizes and P(n)1 and P(n)2 are the search directions which
are determined from:
P(n)1 =
parenleftbigg?J
?p
parenrightbigg(n)
+?(n)1 P(n?1)1 (3.53)
P(n)2 =
parenleftbigg ?J
?Bi
parenrightbigg(n)
+?(n)2 P(n?1)2 . (3.54)
The conjugate coefficients ?(n)1 and ?(n)2 are given by:
?(n)1 =
bracketleftbiggparenleftBig
?J
?p
parenrightBig(n)bracketrightbigg2
bracketleftbiggparenleftBig
?J
?p
parenrightBig(n?1)bracketrightbigg2 with ?
(0)
1 = 0 (3.55)
?(n)2 =
integraltextLe
0
bracketleftbiggparenleftBig
?J
?Bi
parenrightBig(n)bracketrightbigg2
dy
integraltextLe
0
bracketleftbiggparenleftBig
?J
?Bi
parenrightBig(n?1)bracketrightbigg2
dy
with ?(0)2 = 0. (3.56)
The iterative process described by Eq.(3.51)?(3.56) requires the step sizes ?(n)1 and ?(n)2
and the gradient of the functional, ?J(n) =
bracketleftbiggparenleftBig
?J
?p
parenrightBig(n)
,
parenleftBig ?J
?Bi
parenrightBig(n)bracketrightbigg
. These parameters will be
computed from the sensitivity and adjoint problems, respectively.
63
3.3.2 Sensitivity problem and calculation of the search step sizes
The problem of finding the boundary configuration of the computational domain is
a shape identification problem. Shape identification problems were previously studied by
Huang and Shih [23], Huang and Chen [21], Huang and Chao [20], and Park and Shin
[46]. In these works either one or multiple surface configurations were recovered using the
conjugate gradient method. In this analysis the solution of the heat transfer is controlled by
the convective boundary condition, taking the Biot number and the aspect ratio as controls.
We need to solve two sensitivity problems corresponding to each control parameter. The
sensitivity problemsare obtained fromthe direct problem, Eq.(3.1)-(3.14), by perturbingthe
aspect ratio and Biot number, one at a time. When p is perturbed by ?p, the temperatures
TA, TB, and TC are perturbed by ?T1A, ?T1B, and ?T1C, respectively. If Bi undergoes a
perturbation Bi + ?Bi, the temperatures TA, TB, and TC are perturbed by ?T2A, ?T2B,
and ?T2C, respectively. Then replacing in the direct problem p by p+?p, TA by TA+?T1A,
TB by TB + ?T1B, and TC by TC + ?T1C and subtracting from the resulting equations the
direct problem, the following sensitivity problem for the sensitivity functions ?T1A, ?T1B,
and ?T1C is obtained, if the second order terms are neglected:
?2?T1A
?x2 +
?2?T1A
?y2 +
?2?T1A
?z2 = 0 (3.57)
?2?T1B
?x2 +
?2?T1B
?y2 +
?2?T1B
?z2 = 0 (3.58)
?2?T1C
?x2 +
?2?T1C
?y2 +
?2?T1C
?z2 = 0 (3.59)
64
subject to the following boundary conditions,
??T1A
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=0 = 0
??T1A
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=L = 0 (3.60)
??T1A
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=0 = 0
??T1A
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=w
1
= 0 (3.61)
??T1B
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=0 = 0
??T1B
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=W = 0 (3.62)
??T1B
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=0 = 0
??T1B
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=W = 0 (3.63)
??T1C
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=0 = 0
?T1C = TC (l + ?l)?TC (l) ? ?TC?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?l
H1?p = ?Bi TC|x=L?l H1?p
??T1C
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=L?l =
?
???
???
?Bi ?TC?x
vextendsinglevextendsingle
vextendsinglex=L?l H1?p if 0 ? y ? Le
0 if Le < y ? W
(3.64)
??T1C
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=0 = 0
??T1C
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=W = 0 (3.65)
??T1C
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=0 = 0
??T1A
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=1+H
3
= 0 (3.66)
and the conditions of temperature and flux continuity across the GaAs slab-silicon substrate
interface,
?T1B
vextendsinglevextendsingle
vextendsinglez=1 = ?T1A
vextendsinglevextendsingle
vextendsinglez=1 if 0 ? y ? w1 (3.67)
65
??T1B
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=1 =
?
??
???
kGaAs
kSi
??T1A
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
if 0 ? y ? w1
0 if w1 < y ? W
(3.68)
and across the silicon substrate-fin interface,
?T1B
vextendsinglevextendsingle
vextendsinglez=H
1
= ?T1C
vextendsinglevextendsingle
vextendsinglez=H
1
if 0 ? x ? L?l (3.69)
??T1B
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=H
1
=
??
???
???
???
???
??
???
???
???
???
???
??T1C
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=H1
if 0 ? x ? L?l
Bi ?T1Bvextendsinglevextendsinglez=H1 if
?
???
???
L?l < x ? L
0 < y ? Le
0 if
?
???
???
L?l < x ? L
Le < y ? L.
(3.70)
The solution of the sensitivity problem corresponding to a variation of the Biot number
?Bi is obtained by solving the following equations:
?2?T2A
?x2 +
?2?T2A
?y2 +
?2?T2A
?z2 = 0 (3.71)
?2?T2B
?x2 +
?2?T2B
?y2 +
?2?T2B
?z2 = 0 (3.72)
?2?T2C
?x2 +
?2?T2C
?y2 +
?2?T2C
?z2 = 0 (3.73)
with the corresponding boundary conditions,
??T2A
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=0 = 0
??T2A
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=L = 0 (3.74)
66
??T2A
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=0 = 0
??T2A
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=w
1
= 0 (3.75)
??T2B
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=0 = 0
??T2B
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=W = 0 (3.76)
??T2B
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=0 = 0
??T2B
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=W = 0 (3.77)
??T2C
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglex=0 = 0
??T2C
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingleL?l =
?
???
???
???
???
Bi ?T2Cvextendsinglevextendsinglex=L?l + ?Bi TC|x=L?l if 0 ? y ? Le
0 if Le < y ? W
(3.78)
??T2C
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=0 = 0
??T2C
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsingley=W = 0 (3.79)
??T2C
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=0 = 0
??T2A
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=1+H
3
= 0 (3.80)
and the conditions of temperature and flux continuity across the GaAs slab-silicon substrate
interface,
?T2B
vextendsinglevextendsingle
vextendsinglez=1 = ?T2A
vextendsinglevextendsingle
vextendsinglez=1 if 0 ? y ? w1 (3.81)
??T2B
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=1 =
??
??
??
kGaAs
kSi
??T2A
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
if 0 ? y ? w1
0 if w1 < y ? W
(3.82)
and across the silicon substrate-fin interface,
?T2B
vextendsinglevextendsingle
vextendsinglez=H
1
= ?T2C
vextendsinglevextendsingle
vextendsinglez=H
1
if 0 ? x ? L?l (3.83)
67
??T2B
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
vextendsinglez=H
1
=
?
???
???
???
???
???
???
???
???
???
???
??T2C
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=H1
if 0 ? x ? L?l
Bi ?T2Bvextendsinglevextendsinglez=H1 + ?Bi TB|z=H1 if
?
???
???
L?l < x ? L
0 < y ? Le
0 if
?
???
???
L?l < x ? L
Le < y ? L.
(3.84)
The functional J(n+1) for the (n+1)th iteration is written by replacing p(n+1) and Bi(n+1)
with the expressions given by Eq.(3.51)-(3.52):
J(n+1) =
integraldisplay 1+H3
1
integraldisplay w1
0
integraldisplay L
0
bracketleftBig
TA
parenleftBig
x,y,z,p(n) ??(n)1 P(n)1 ,Bi(n) ??(n)2 P(n)2
parenrightBigbracketrightBig2
dxdydz
+?2
integraldisplay Le
0
parenleftBig
Bi(n) ??(n)2 P(n)2
parenrightBig2
dy. (3.85)
By expanding in Taylor series TA, p, and Bi the following expression is obtained for the
functional J at the (n+ 1)th iteration:
J(n+1) =
integraldisplay 1+H3
1
integraldisplay w1
0
integraldisplay L
0
bracketleftBig
TA
parenleftBig
x,y,z,p(n),Bi(n)
parenrightBig
??(n)1 ?T1A
parenleftBig
P(n)1
parenrightBig
? ?(n)2 ?T2A
parenleftBig
P(n)2
parenrightBigbracketrightBig2
dxdydz + ?2
integraldisplay Le
0
parenleftBig
Bi(n) ??(n)2 P(n)2
parenrightBig2
dy. (3.86)
The search step sizes ?1 and ?2 are found from the line search in each direction at each
iteration by taking partial derivatives of Eq.(3.86) with respect to ?1 and ?2 and equating
the obtained expressions to zero. The following system of equations is obtained:
?(n)1 =
integraltext1+H3
1
integraltextw1
0
integraltextL
0 2?T
1AparenleftbigTA ??2?T2Aparenrightbigdxdydz
integraltext1+H3
1
integraltextw1
0
integraltextL
0 2
parenleftbig?T1
A
parenrightbig2 dxdydz (3.87)
68
?(n)2 =
integraltext1+H3
1
integraltextw1
0
integraltextL
0 2?T
2AparenleftbigTA ??1?T1Aparenrightbigdxdydz + ?
2
integraltextLe
0 Bi
(n)P(n)2 dy
integraltext1+H3
1
integraltextw1
0
integraltextL
0 2
parenleftbig?T2
A
parenrightbig2 dxdydz + ?
2
integraltextLe
0
parenleftBig
P(n)2
parenrightBig2
dy
. (3.88)
3.3.3 Adjoint problems
The Lagrangian associated with the optimization problem defined by Eq.(3.50), the
constraints (3.1)-(3.14), and the control parameter p is:
J|p =
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
T2Adxdydz + ?2
integraldisplay Le
0
Bi2dy
+
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
?1
parenleftBigg
?2TA
?x2 +
?2TA
?y2 +
?2TA
?z2
parenrightBigg
dxdydz
+
integraldisplay 1
H1
integraldisplay W
0
integraldisplay L
0
?2
parenleftBigg
?2TB
?x2 +
?2TB
?y2 +
?2TB
?z2
parenrightBigg
dxdydz
+
integraldisplay H1
0
integraldisplay W
0
integraldisplay L?pH1
0
?3
parenleftBigg
?2TC
?x2 +
?2TC
?y2 +
?2TC
?z2
parenrightBigg
dxdydz (3.89)
where ?1, ?2, and ?3 are Lagrange multipliers. The differential of Lagrangian corresponding
to ?p is:
J|p+?p ? J|p =
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
2TA?T1Adxdydz
+
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
?1
parenleftBigg
?2?T1A
?x2 +
?2?T1A
?y2 +
?2?T1A
?z2
parenrightBigg
dxdydz
+
integraldisplay 1
H1
integraldisplay W
0
integraldisplay L
0
?2
parenleftBigg
?2?T1B
?x2 +
?2?T1B
?y2 +
?2?T1B
?z2
parenrightBigg
dxdydz
+
integraldisplay H1
0
integraldisplay W
0
integraldisplay L?H1p
0
?3
parenleftBigg
?2?T1C
?x2 +
?2?T1C
?y2 +
?2?T1C
?z2
parenrightBigg
dxdydz. (3.90)
69
Integration of the above expression leads to the following adjoint problem after boundary
conditions from the sensitivity problem (3.57)-(3.70) are applied:
?2?1
?x2 +
?2?1
?y2 +
?2?1
?z2 + 2TA = 0 (3.91)
?2?2
?x2 +
?2?2
?y2 +
?2?2
?z2 = 0 (3.92)
?2?3
?x2 +
?2?3
?y2 +
?2?3
?z2 = 0 (3.93)
subject to the following boundary conditions,
??1
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??1?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L
= 0 (3.94)
??1
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??1?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=w1
= 0 (3.95)
??2
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??2?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?l
= 0 (3.96)
??2
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??2?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=W
= 0 (3.97)
??3
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??3?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?l
= 0 (3.98)
??3
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??3?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=W
= 0 (3.99)
??3
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=0
= 0 ??1?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1+H3
= 0 (3.100)
70
and the conditions of temperature and flux continuity across the GaAs slab-silicon substrate
interface,
?2|z=1 = kGaAsk
Si
?1|z=1 if 0 ? y ? w1 (3.101)
??2
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
=
??
??
???
??1
?z
vextendsinglevextendsingle
vextendsinglez=1 if 0 ? y ? w1
0 if w1 < y ? W
(3.102)
and across the silicon substrate-fin interface,
?2|z=H1 = ?3|z=H1 if 0 ? x ? L?l (3.103)
??2
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=H1
=
?
???
???
???
???
??
???
???
???
???
??
??2
?z
vextendsinglevextendsingle
vextendsinglez=H
1
if 0 ? x ? L?l
Bi?2|z=H1 if
??
??
???
L?l < x ? L
0 ? y ? Le
0 if
?
???
???
L?l < x ? L
Le < y ? W.
(3.104)
After identifying the adjoint problem, the following integral term is left:
?J|p = J|p+?p ? J|p
=
integraldisplay H1
0
integraldisplay Le
0
Bi?TC?x ?3
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?pH1
?pdydz
=
integraldisplay H1
0
integraldisplay Le
0
?J
?p?pdydz (3.105)
71
and the following expression for ?J?p is obtained,
?J
?p = Bi
?TC
?x ?3
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?pH1
. (3.106)
The adjoint problem corresponding to Bi is obtained in a similar way. The associated
Lagrangian in this case is:
J|Bi =
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
T2Adxdydz + ?2
integraldisplay Le
0
Bi2dy
+
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
?1
parenleftBigg
?2TA
?x2 +
?2TA
?y2 +
?2TA
?z2
parenrightBigg
dxdydz
+
integraldisplay 1
H1
integraldisplay W
0
integraldisplay L
0
?2
parenleftBigg
?2TB
?x2 +
?2TB
?y2 +
?2TB
?z2
parenrightBigg
dxdydz
+
integraldisplay H1
0
integraldisplay W
0
integraldisplay L?pH1
0
?3
parenleftBigg
?2TC
?x2 +
?2TC
?y2 +
?2TC
?z2
parenrightBigg
dxdydz (3.107)
and the differential of Lagrangian corresponding to ?Bi is:
J|Bi+?Bi ? J|Bi =
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
2TA?T2Adxdydz +?
integraldisplay Le
0
Bi?Bidy
+
integraldisplay 1+H3
1
integraldisplay W
0
integraldisplay L
0
?1
parenleftBigg
?2?T2A
?x2 +
?2?T2A
?y2 +
?2?T2A
?z2
parenrightBigg
dxdydz
+
integraldisplay 1
H1
integraldisplay W
0
integraldisplay L
0
?2
parenleftBigg
?2?T2B
?x2 +
?2?T2B
?y2 +
?2?T2B
?z2
parenrightBigg
dxdydz
+
integraldisplay H1
0
integraldisplay W
0
integraldisplay L?pH1
0
?3
parenleftBigg
?2?T2C
?x2 +
?2?T2C
?y2 +
?2?T2C
?z2
parenrightBigg
dxdydz.
(3.108)
72
The following adjoint problem is obtained after performing the integration by parts in
Eq.(3.108):
?2?1
?x2 +
?2?1
?y2 +
?2?1
?z2 + 2TA = 0 (3.109)
?2?2
?x2 +
?2?2
?y2 +
?2?2
?z2 = 0 (3.110)
?2?3
?x2 +
?2?3
?y2 +
?2?3
?z2 = 0 (3.111)
subject to the following boundary conditions,
??1
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??1?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L
= 0 (3.112)
??1
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??1?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=w1
= 0 (3.113)
??2
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??2?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L
= 0 (3.114)
??2
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??2?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=W
= 0 (3.115)
??3
?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=0
= 0 ??3?x
vextendsinglevextendsingle
vextendsinglevextendsingle
x=L?l
=
??
??
???
?Bi ?3|x=L?l if 0 ? y ? Le
0 if Le < y ? W
(3.116)
??3
?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=0
= 0 ??3?y
vextendsinglevextendsingle
vextendsinglevextendsingle
y=W
= 0 (3.117)
??3
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=0
= 0 ??1?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1+H3
= 0 (3.118)
and the conditions of temperature and flux continuity across the GaAs slab-silicon substrate
interface,
?1|z=1 = kGaAsk
Si
?2|z=1 if 0 ? y ? w1 (3.119)
73
??2
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=1
=
?
??
???
??1
?z
vextendsinglevextendsingle
vextendsinglez=1 if 0 ? y ? w1
0 if w1 < y ? W
(3.120)
and across the silicon substrate-fin interface,
?2|z=H1 = ?3|z=H1 if 0 ? x ? L?l (3.121)
??2
?z
vextendsinglevextendsingle
vextendsinglevextendsingle
z=H1
=
?
???
???
???
???
??
???
???
???
???
??
??3
?z
vextendsinglevextendsingle
vextendsinglez=H
1
if 0 ? x ? L?l
Bi ?2|z=H1 if
?
???
???
L?l < x ? L
0 < y ? Le
0 if
?
???
???
L?l < x ? L
Le < y ? L.
(3.122)
After identifying the adjoint problem, the following integral term is left:
?J|Bi = J|Bi+?Bi ? J|Bi
=
integraldisplay Le
0
parenleftBigg
?Bi?
integraldisplay H1
0
?3TC|x=L?pH1 dz ?
integraldisplay L
L?pH1
?2TB|z=H1 dx
parenrightBigg
?Bidy
(3.123)
and the following expression for ?J?Bi is obtained,
?J
?Bi = ?Bi?
integraldisplay H1
0
?3TC|x=L?pH1 dz ?
integraldisplay L
L?pH1
?2TB|z=H1 dx. (3.124)
74
3.3.4 Computational algorithm
The computational procedure for the solution of this problem may be summarized as
follows:
step 1 Solve the direct problem given by Eq.(3.1)-(3.14) to obtain TA, TB and TC for an
initial guess of the controls p and Bi;
step 2 Solve the adjoint problems given by Eq.(3.91)-(3.104) and Eq.(3.109)-(3.122), respec-
tively to obtain ?1 ?2, ?3, ?1, ?2, and ?3;
step 3 Compute the gradients of the functional ?J?p and ?J?Bi from Eq.(3.106) and Eq.(3.124),
respectively;
step 4 Compute the conjugate coefficients ?(n)1 and ?(n)2 from Eq.(3.55) and Eq.(3.56), and
the search directions P(n)1 and P(n)2 from Eq.(3.53) and Eq.(3.54);
step 5 Set ?p = P(n)1 , ?Bi = P(n)2 and solve the sensitivity problems given by Eq.(3.57)-
(3.70) and Eq.(3.71)-(3.84) to obtain ?T1A, ?T1B, ?T1C and ?T2A, ?T2B, ?T2C, respec-
tively;
step 6 Compute the search step sizes ?(n)1 and ?(n)2 from Eq.(3.87) and Eq.(3.88);
step 7 Compute the new estimates for p and Bi from Eq.(3.51) and Eq.(3.52) and go back
to step 1 until a stopping criterion is achieved. The stopping criterion was chosen as
vextenddoublevextenddouble
vextenddoublep(n+1) ?p(n)
vextenddoublevextenddouble
vextenddouble? ?1,
vextenddoublevextenddouble
vextenddoubleBi(n+1) ?Bi(n)
vextenddoublevextenddouble
vextenddouble? ?2, where ?1 and ?2 are prescribed tolerances
which were set to 10?6 and 10?7, respectively.
75
3.3.5 Implementation and numerical procedure
The direct, adjoint, and sensitivity problems are solved using a centered finite difference
scheme. The following geometrical and thermophysical parameters are given:
? GaAs dimensions: Lt ?w1 ?H3 = 9mm?3.6mm?0.075mm;
? Si wafer dimensions: Lt ?W ?H2 = 9mm?9mm?0.3mm;
? Coefficient of thermal conductivity for GaAs: kGaAs = 0.1W/mmK;
? Coefficient of thermal conductivity for Si: kSi = 0.1421W/mmK;
? Input heat flux: q = 12W/mm2.
A Matlab code is implemented to solve the three coupled problems using a mesh of 15?11?4
for the GaAs slab and 15?25?12 for the Si substrate.
3.3.6 Results and discussions
The optimization of the MHP sink is performed for 4 different channel heights, H1 =
75?m,100?m,125?m,150?m. The number of channels is also varied to study the sensitivity
of the heat sink performance on the number of channels. The optimum channel to fin width
ratio is 1, this finding being in agreement with the results obtained by previous investigators
[68], [32], [53]. The obtained aspect ratios result in reasonable fin widths and arrays of
fins that can be manufactured. Distributions of optimum Bi are found for each case study,
considering five values of the weighting coefficient, ? = 0.1,0.5,1,5,50. It is noticed that the
optimum Biot number distribution decreases with increasing channel depth and increases
with decreasing number of channels (see Figures 4.5-4.8). As in the planar heat sink case
studied in the previous chapter, Bi decreases with the increase of cost ?.
76
The controlled and uncontrolled temperature distributions on the top surface of the Si
substrate through the symmetry plane of the fin are compared for different values of the
weighting coefficient ? (see Figures 4.9-4.12). The uncontrolled temperature is obtained
using a constant Biot number, BI along the boundary, that gives the same cost as the
optimal Bi. It can be noticed that the controlled temperatures display flatter profiles,
resulting in smaller temperature gradients. Underneath the heat source, the controlled
temperature is lower than the uncontrolled temperature. The difference between the two
temperatures corresponding to the hottest point is about 4K ? 10K. This difference is
much smaller than in the planar spreader case, suggesting that the heat transfer coefficient
of a finned structure approaches by design an optimum distribution. Since the differences
between the controlled and uncontrolled temperatures are small, the constant Biot number
BI is used to assess the effect of the numberof channels on the performanceof the MHP sink.
Thestructure withN = 60 channels is taken as reference and variation percentages of BI are
computed when the number of channels is changed to N = 58, 56 and 30, respectively (see
Table 4.1). For a given number of channels no significant variations are noticed for different
values of ?. As can be observed from Figures 4.9-4.12, for a given channel depth H1 and
weighting coefficient ? the same temperature distributions (controlled and uncontrolled) are
obtained for values of the Biot number that increase with decreasing number of channels.
For the same number of channels lower temperatures are obtained for larger values
of the channel depth. By increasing the channel depth the surface exposed to convective
cooling is expanded.
Comparing the thermal performance of the planar spreader studied in the previous
chapter and the finned heat sinks, it can be concluded that the same temperatures can be
77
attained for smaller Bi using the finned heat sinks. This is due to an increase in area for
heat dissipation.
3.4 Liquid charge corresponding to maximum heat transport capacity of the
pipe
In this section we will compute the liquid charges corresponding to the maximum heat
transport capacity of the MHPs having the optimum geometry found in the previous section.
It is expected that the optimum Bi found in the previous section leads to an amount of
heat removed from the source greater than the heat transport capacity of the MHP. The
maximum heat transport capacity of the pipe is found iteratively, decreasing the optimum
Bi by a quantity ?Bi(i) = i?10?7 ?Bi, i = 0..107 .
Q(i)inMHP = qin (H1 +l)
integraldisplay Le
0
parenleftBig
Bi??Bi(i)
parenrightBig
Twdy. (3.125)
In the above equation qin is the heat flux entering the heat sink and Tw is the channel wall
temperature. The maximum heat applied to the MHP for which the pressure gradients
dPl
dy ,
dPv
dy do not change signs along the evaporator and adiabatic sections, and the meniscus
radius at the end of the adiabatic section does not exceed the rmax value (see Eq.(3.27)) is
the maximum heat transport capacity of the MHP. A new optimum Bi is computed taking
the maximum heat transport capacity of the pipe as the input heat for the heat sink. The
liquid charge is found by adding the amount of liquid in the evaporator, adiabatic section,
and condenser and the amount of liquid that vaporizes:
78
mliquid =
parenleftBigg
2
integraldisplay Le+La
0
Aldy +2lH1 (W ?Le ?La)
parenrightBigg
?l
+
parenleftBigg
2lH1 (Le +La)?2
integraldisplay Le+La
0
Aldy
parenrightBigg
?v. (3.126)
3.4.1 Computational procedure
The algorithm used in computing the optimum liquid charge may be summarized as
follows:
step 1 For a given r|y=0 = r0 set the heat entering the MHP to Q0inMHP (i=0);
step 2 Solve the transport Eq.(3.44)-(3.48) for the MHP to obtain the distributions for Pl,
Pv, Ul, Uv, and r;
step 3 Check the conditions for the pressure gradients and the radius of curvature:
dPl
dy < 0,
dPv
dy > 0, r ? rcrit y ? Le +La; (3.127)
step 4 If all of the above conditions are satisfied, then the heat transport capacity of the
MHP is Qmax = QiinMHP, otherwise increase ?Bi (set the input heat to Q(i+1)inMHP)
and repeat step 2-step 4 until the above conditions are satisfied simultaneously;
step 5 Take Qmax as the input heat for the MHP heat sink and compute the optimum
Bi using a conjugate gradient algorithm, similar to those presented in the previous
sections;
step 6 Compute the liquid charge using Eq.(3.126).
79
3.4.2 Results and discussions
In this analysis water was considered as working fluid. The following geometrical and
thermophysical parameters are used:
? MHP height: H1 = 75?m,100?m,125?m,150?m;
? MHP width: l = pH1, where p is the optimum aspect ratio found in the previous
section;
? MHP length: W = 9mm;
? Evaporator length: Le = 1.8mm;
? Adiabatic section length: La = 3mm;
? Liquid density: ?l = 983.284kg/m3;
? Vapor density: ?v = 0.130366kg/m3
? Liquid viscosity: ?l = 466.4 ?10?6kg/ms;
? Vapor viscosity: ?v = 10.93 ?10?6kg/ms;
? Surface tension: ? = 0.06624N/m;
? Contact angle: ? = ?/6;
? Tilt angle: ?1 = pi/18
? Latent heat of vaporization: hfg = 2358.48kJ/kg.
80
The properties for water are taken at 55?C.
The transport equations (3.44)-(3.48) for the MHP are solved using a Runge-Kutta
scheme. Distributions of the meniscus radius over the axial direction for the evaporator
and adiabatic sections are shown in Figure 4.13. These distributions correspond to the
maximum heat transport capacity of MHP sinks with N = 60 channels and different channel
depths. It is noticed that the meniscus radius increases with increasing channel depth.
Figures 4.14-4.17 show the velocity distributions of the liquid and vapor phases computed
for the maximum heat transport capacity of the above mentioned heat sink designs. Higher
velocities are achieved for greater channel depths. The increase of radius of curvature in
the axial direction causes an increase in liquid cross section area in this direction and a
decrease in vapor area. Consequently the vapor velocity increases in y direction while the
liquid velocity decreases. The vapor velocity is about two orders of magnitude greater than
the liquid velocity. This is due to the high liquid/vapor density ratio. Pressure distributions
for the liquid and vapor phases are presented in Figure 4.18. The vapor pressure drop is
about one order of magnitude less than that of the liquid phase.
Optimal distributions of Bi corresponding to the maximum heat transport capacity
of the MHP are presented in Figures 4.19-4.22 for heat sinks with different number of
channels and different channel depths. The new optimum Bi represents only a fraction
of the optimum Bi computed in the previous section (?ideal? Bi). For a given number of
channels the optimum Bi increases with increasing channel depth H1. Values of the liquid
charge corresponding to the maximum heat transport capacity of the pipe are plotted in
Figure 4.25. For a given number of channels the optimum liquid charge increases with
increasing channel depth. Small increases in the liquid charge are noticed for increasing
81
weighting coefficient ?. The optimum liquid charge increases with decreasing number of
channels.
The maximum heat transport capacity is expressed as a fraction qr of the input heat.
Variations of the MHP heat transport capacity for different design parameters are shown
in Figure 4.27. For a given number of channels, the heat transport capacity increases with
increasing channel depth and weighting coefficient ?.
The sensitivity of the MHP performance on the value of the meniscus radius at y = 0 is
also studied. Four values of r0 were considered, 7.7?m,8?m,8.5?m and 9.5?m. The results
obtained for the optimum liquid charge (see Figure 4.23) and maximum heat transport
capacity (see Figure 4.24) do not display significant variations when r0 is varied within
25%.
3.5 Conclusions
In this work the performance of several configurations of rectangular cross section
micro heat pipe sinks using water as working fluid was investigated. The geometry and
the convective heat transfer coefficient were optimized. The critical path to an optimum
design started with the analysis of a planar heat spreader for which an optimum convective
heat transfer coefficient corresponding to a maximum amount of heat removed from a heat
source was obtained. An optimal control technique was employed to solve this problem,
using Biot number as control. The existence and uniqueness of a solution and of an optimal
control were proven. Then the analysis was extended to finned heat sinks for which an
optimization of the convective heat transfer coefficient and the aspect ratio was performed.
82
The last part of this study focused on obtaining the amount of liquid charge corresponding
to the maximum heat transport capacity of the micro heat pipes.
For the planar spreader, distributions of the Biot number corresponding to different
values of the weighting coefficient ? were presented, showing that Bi decreases with in-
creasing ?. The greater the value assigned to ?, the greater the cost of employing a control,
and the smaller the contribution of convective cooling. Controlled and uncontrolled tem-
perature distributions were compared. The controlled temperatures displayed lower and
flatter profiles than the uncontrolled temperatures. The difference between the two profiles
increased with increasing penalty parameter ?.
A similar optimal control technique was used to optimize the geometry and the con-
vective coefficient of finned heat sinks. The results show that the optimum aspect ratio
corresponds to a channel to fin width ratio equal to 1. The distributions of the optimum
Biot number are lower than in the case of a planar spreader, due to an increase in heat
transfer area. An analysis of heat sink configurations with different channel depths, i.e.
75?m, 100?m, 125?m, and 150?m, showed that the performance of the heat sink increases
with increasing channel depth. The performance sensitivity on the number of channels was
studied. The sink with 60 channels was taken as reference and computations of the varia-
tion percentage of Biot number were performed when the number of channels was changed
to 58, 56, and 30. When the number of channels was reduced to half, the Biot number
was increased by 35%-43%, suggesting that the performance of a heat sink increases with
increasing channel density.
The maximum heat transport capacities of MHP sinks with an optimum geometry were
computed, and found to be about 5%?15% of the original input heat flux of 12W/mm2, i.e.
83
between 0.6W/mm2 and 1.8W/mm2. The heat transport capacity increases with increasing
channel density and increasing aspect ratio. The liquid charge found was about 9%?11%
(see Figure 4.26) of the pipe volume.
The algorithm used in computing the liquid charge can be applied to any working
fluid. The only change that has to be made in the simulation consists in replacing the
thermophysical properties of water with the ones corresponding to the new working fluid.
For applications where the removal of higher heat fluxes is necessary, different working fluids
have to be used, such as mercury [44].
This study can provide guidance and a methodology for design of efficient heat sinks
using micro heat pipes. The design procedure can be followed as given:
? Choose the depth of the micro-channels. Typical values are in the range of 100?m,
but can vary from 50mm up to several hundred microns. The determining factor lies
in the thickness of the spreader and the location of heat generation. In general, the
channels should be as close to the heat source as possible.
? For the given dimensions of the heat sink, length by width (with the width being
defined as the dimension transverse of the direction of the heat pipe channel axial
dimension) prescribe the number of channels and find the channel width using the
recommended pitch ratio (channel to fin width ratio) of 1 as described in section 3.3,
e.g. for a 10mm (width) by 10mm (length) heat spreader with a 0.5mm border on
all edges and 100?m deep channels, if the recommended number of channels is 60
across the width, then the channel width is 75?m. The number of channels is selected
considering the amount of heat that has to be dissipated. For example, the above
mentioned configuration is able to handle a heat flux of 0.9W/mm2 (see Figure 4.27).
84
? The corresponding liquid charge can be found using the algorithm described in section
3.4.1. This value is somewhere between 9 and 11%.
As designers we look for a heat sink configuration able to spread the largest amount of heat
for a particular application. From the study presented earlier the best performance was
achieved by a heat sink with 60 channels using a channel depth of 150?m, which was able
to handle a heat flux of 1.8W/mm2 for a 11% volume charge.
85
Chapter 4
Appendix
Table 4.1: Variation percentage of BI for different number of channels
H1 = 75?m H1 = 100?m H1 = 125?m H1 = 150?m
N = 58 1.8124% 1.9558% 2.1194% 2.2216%
N = 56 3.8110% 3.9912% 4.1710% 4.3048%
N = 30 34.9441% 39.3297% 40.7426% 43.2075%
86
0 1 2 3 4 5 6 7 8 90
50
100
150
200
250
300
x [mm]
T?T
? [K]
Centerline temperature distribution, ?=0.1
controlled solution
uncontrolled solution
0 1 2 3 4 5 6 7 8 90
50
100
150
200
250
300
x [mm]
T?T
? [K]
Centerline temperature distribution, ?=0.5
controlled solution
uncontrolled solution
0 1 2 3 4 5 6 7 8 90
50
100
150
200
250
300
x [mm]
T?T
? [K]
Centerline temperature distribution, ?=1
controlled solution
uncontrolled solution
0 1 2 3 4 5 6 7 8 90
50
100
150
200
250
300
x [mm]
T?T
? [K]
Centerline temperature distribution, ?=5
controlled solution
uncontrolled solution
0 1 2 3 4 5 6 7 8 90
50
100
150
200
250
300
x [mm]
T?T
? [K]
Centerline temperature distribution, ?=50
controlled solution
uncontrolled solution
Figure 4.1: Graphs of controlled and uncontrolled solutions for different values of ?.
87
0
3
6
9
0
3
6
9
0
0.5
1
1.5
2
2.5
3
x [mm]
Control function distribution, ?=0.1
y [mm]
Bi
0.5
1
1.5
2
2.5
0
3
6
9
0
3
6
9
0
0.5
1
1.5
2
2.5
3
x [mm]
Control function distribution, ?=0.5
y [mm]
Bi
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
0
3
6
9
0
3
6
9
0
0.5
1
1.5
2
2.5
3
x [mm]
Control function distribution, ?=1
y [mm]
Bi
0.2
0.4
0.6
0.8
1
1.2
1.4
0
3
6
9
0
3
6
9
0
0.5
1
1.5
2
2.5
3
x [mm]
Control function distribution, ?=5
y [mm]
Bi
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0
3
6
9
0
3
6
9
0
0.5
1
1.5
2
2.5
3
x [mm]
Control function distribution, ?=50
y [mm]
Bi
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Figure 4.2: Distributions of Bi for different values of ?.
88
0
2
4
6
8
0
2
4
6
8
0
50
100
150
x [mm]
Temperature distribution at the interface, ?=0.1
y [mm]
T?T
? [K]
5
10
15
20
25
30
35
40
45
50
0
3
6
9
0
3
6
9
0
50
100
150
x [mm]
Temperature distribution at the interface, ?=0.5
y [mm]
T?T
? [K]
10
20
30
40
50
60
0
3
6
9
0
3
6
9
0
50
100
150
x [mm]
Temperature distribution at the interface, ?=1
y [mm]
T?T
? [K]
10
20
30
40
50
60
0
3
6
9
0
3
6
9
0
50
100
150
x [mm]
Temperature distribution at the interface, ?=5
y [mm]
T?T
? [K]
10
20
30
40
50
60
70
80
0
3
6
9
0
3
6
9
0
50
100
150
x [mm]
Temperature distribution at the interface, ?=50
y [mm]
T?T
? [K]
20
40
60
80
100
120
Figure 4.3: Distributions of T ?T? at the chip-spreader interface for different values of ?.
89
0 1 2 3 4 5 6 7 8 90
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
x [mm]
z [mm]
Temperature distribution through the central cross section, ?=0.1
10
20
30
40
50
60
0 1 2 3 4 5 6 7 8 90
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
x [mm]
z [mm]
Temperature distribution through the central cross section, ?=0.5
10
20
30
40
50
60
0 1 2 3 4 5 6 7 8 90
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
x [mm]
z [mm]
Temperature distribution through the central cross section, ?=1
58
60
62
64
66
68
70
72
74
0 1 2 3 4 5 6 7 8 90
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
x [mm]
z [mm]
Temperature distribution through the central cross section, ?=5
72
74
76
78
80
82
84
86
88
90
92
0 1 2 3 4 5 6 7 8 90
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
x [mm]
z [mm]
Temperature distribution through the central cross section, ?=50
105
110
115
120
125
Figure 4.4: Distributions of T ?T? through a vertical symmetry plane for different values
of ?.
90
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=60, H1=75?m, pcrit=0.9333
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=60, H1=100?m, pcrit=0.7
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=60, H1=125?m, pcrit=0.56
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=60, H1=150?m, pcrit=0.4667
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
Figure 4.5: Distributions of optimal Bi corresponding to N = 60 and different channel
depths.
91
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=58, H1=75?m, pcrit=0.9678
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=58, H1=100?m, pcrit=0.7259
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=58, H1=125?m, pcrit=0.5807
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=58, H1=150?m, pcrit=0.4839
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
Figure 4.6: Distributions of optimal Bi corresponding to N = 58 and different channel
depths.
92
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=56, H1=75?m, pcrit=1.0048
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=56, H1=100?m, pcrit=0.7536
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=56, H1=125?m, pcrit=0.6029
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=56, H1=150?m, pcrit=0.5024
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
Figure 4.7: Distributions of optimal Bi corresponding to N = 56 and different channel
depths.
93
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=30, H1=75?m, pcrit=1.9333
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=30, H1=100?m, pcrit=1.45
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=30, H1=125?m, pcrit=1.16
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
N=30, H1=150?m, pcrit=0.9667
y [mm]
Bi
?=0.1
?=0.5
?=1
?=5
?=50
Figure 4.8: Distributions of optimal Bi corresponding to N = 30 and different channel
depths.
94
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=60, H1=75?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=60, H1=100?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=60, H1=125?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=60, H1=150?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
Figure 4.9: Temperature distributions T ?T? corresponding to N = 60.
95
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=58, H1=75?m ?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=58, H1=100?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=58, H1=125?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=58, H1=150?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
Figure 4.10: Temperature distributions T ?T? corresponding to N = 58.
96
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=56, H1=75?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=56, H1=100?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=56, H1=125?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=56, H1=150?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
Figure 4.11: Temperature distributions T ?T? corresponding to N = 56.
97
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=30, H1=75?m ?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=30, H1=100?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=30, H1=125?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
10
20
30
40
50
60
70
80
90
100
110
y [mm]
T?T
? [K]
N=30, H1=150?m
?=0.1, controlled
?=0.1, uncontrolled
?=0.5, controlled
?=0.5, uncontrolled
?=1, controlled
?=1, uncontrolled
?=5, controlled
?=5, uncontrolled
?=50, controlled
?=50, uncontrolled
Figure 4.12: Temperature distributions T ?T? corresponding to N = 30.
98
0 0.5 1 1.5 2 2.5 3 3.50
5
10
15
20
25
30
35
40
45
y[mm]
r[?
m]
H1=75?m
?=0.1
?=0.5
?=1
?=5
?=50
0 0.5 1 1.5 2 2.5 3 3.50
5
10
15
20
25
30
35
40
45
y[mm]
r[?
m]
H1=100?m
?=0.1
?=0.5
?=1
?=5
?=50
0 0.5 1 1.5 2 2.5 3 3.50
5
10
15
20
25
30
35
40
45
y[mm]
r[?
m]
H1=125?m
?=0.1
?=0.5
?=1
?=5
?=50
0 0.5 1 1.5 2 2.5 3 3.50
5
10
15
20
25
30
35
40
45
y[mm]
r[?
m]
H1=150?m
?=0.1
?=0.5
?=1
?=5
?=50
Figure 4.13: Distribution of the meniscus radius for different values of ? and channel depths
(N = 60).
99
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=75?m, ?=0.1
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=75?m, ?=0.5
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=75?m, ?=1
3.80
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=75?m, ?=5
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=75?m, ?=50
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
Figure 4.14: Distributions of the liquid and vapor velocities for H1 = 75?m and different
values of ? (N = 60).
100
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=100?m, ?=0.1
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=100?m, ?=0.5
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=100?m, ?=1
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=100?m, ?=5
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=100?m, ?=50
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
Figure 4.15: Distributions of the liquid and vapor velocities for H1 = 100?m and different
values of ? (N = 60).
101
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=125?m, ?=0.1
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=125?m, ?=0.5
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=125?m, ?=1
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=125?m, ?=5
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=125?m, ?=50
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
Figure 4.16: Distributions of the liquid and vapor velocities for H1 = 125?m and different
values of ? (N = 60).
102
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=150?m, ?=0.1
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=150?m, ?=0.5
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=150?m, ?=1
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=150?m, ?=5
3.80
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
0 0.5 1 1.5 2 2.5 3 3.5?1
?0.9
?0.8
?0.7
?0.6
?0.5
?0.4
?0.3
?0.2
?0.1
0
U l[m/s]
y[mm]
H1=150?m, ?=50
0
10
20
30
40
50
60
70
80
90
100
U v
[m/s]
Figure 4.17: Distributions of the liquid and vapor velocities for H1 = 150?m and different
values of ? (N = 60).
103
0 0.5 1 1.5 2 2.5 3 3.51.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2 x 10
4
y[mm]
P l,P
v [Pa]
H1=75?m
Pv, ?=0.1
Pl, ?=0.1
Pv, ?=0.5
Pl, ?=0.5
Pv, ?=1
Pl, ?=1
Pv, ?=5
Pl, ?=5
Pv, ?=50
Pl, ?=50
0 0.5 1 1.5 2 2.5 3 3.51.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2 x 10
4
y[mm]
P l,P
v [Pa]
H1=100?m
Pv, ?=0.1
Pl, ?=0.1
Pv, ?=0.5
Pl, ?=0.5
Pv, ?=1
Pl, ?=1
Pv, ?=5
Pl, ?=5
Pv, ?=50
Pl, ?=50
0 0.5 1 1.5 2 2.5 3 3.51.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2 x 10
4
y[mm]
P l,P
v [Pa]
H1=125?m
Pv, ?=0.1
Pl, ?=0.1
Pv, ?=0.5
Pl, ?=0.5
Pv, ?=1
Pl, ?=1
Pv, ?=5
Pl, ?=5
Pv, ?=50
Pl, ?=50
0 0.5 1 1.5 2 2.5 3 3.51.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
2 x 10
4
y[mm]
P l,P
v [Pa]
H1=150?m
Pv, ?=0.1
Pl, ?=0.1
Pv, ?=0.5
Pl, ?=0.5
Pv, ?=1
Pl, ?=1
Pv, ?=5
Pl, ?=5
Pv, ?=50
Pl, ?=50
Figure 4.18: Distributions of the liquid and vapor pressures for different values of ? and
channel depths (N = 60).
104
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=75?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=100?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=125?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=150?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
Figure 4.19: Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding to
the maximum heat transport capacity of the pipe, for N = 60 and different channel depths.
105
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=75?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=100?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=125?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=150?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
Figure 4.20: Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding to
the maximum heat transport capacity of the pipe, for N = 58 and different channel depths.
106
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=75?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=100?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=125?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=150?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
Figure 4.21: Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding to
the maximum heat transport capacity of the pipe, for N = 56 and different channel depths.
107
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=75?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=100?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=125?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.80
0.2
0.4
0.6
0.8
1
1.2
y[mm]
Bi
H1=150?m
?=0.1
?=0.5
?=1
?=5
?=50
optimum Bi
MHP capacity
Figure 4.22: Distributions of the ?ideal? optimum Bi and the optimum Bi corresponding to
the maximum heat transport capacity of the pipe, for N = 30 and different channel depths.
108
70 80 90 100 110 120 130 140 1505
6
7
8
9
10
11
?=0.1
H1[?m]
V liquid
[nl]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1505
6
7
8
9
10
11
?=0.5
H1[?m]
V liquid
[nl]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1505
6
7
8
9
10
11
?=1
H1[?m]
V liquid
[nl]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1505
6
7
8
9
10
11
?=5
H1[?m]
V liquid
[nl]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1505
6
7
8
9
10
11
12
?=50
H1[?m]
V liquid
[nl]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
Figure 4.23: Variations of the liquid charge with r0 for different channel depths (N = 60).
109
70 80 90 100 110 120 130 140 1504
5
6
7
8
9
10
11
12
13
14
?=0.1
H1[?m]
q r[%]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1504
5
6
7
8
9
10
11
12
13
14
?=0.5
H1[?m]
q r[%]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1504
5
6
7
8
9
10
11
12
13
14
?=1
H1[?m]
q r[%]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1504
5
6
7
8
9
10
11
12
13
14
?=5
H1[?m]
q r[%]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
70 80 90 100 110 120 130 140 1504
6
8
10
12
14
16
?=50
H1[?m]
q r[%]
r0=7.7?m r0=8?m r0=8.5?m r0=9.5?m
Figure 4.24: Variations of MHP heat transport capacity with r0 for different channel depths
(N = 60).
110
10?1 100 101 102
4
6
8
10
12
14
16
18
20
?
V liquid
[nl]
N=60
H1=75?m H1=100?m H1=125?m H1=150?m
10?1 100 101 102
4
6
8
10
12
14
16
18
20
?
V liquid
[nl]
N=58
H1=75?m H1=100?m H1=125?m H1=150?m
10?1 100 101 102
4
6
8
10
12
14
16
18
20
?
V liquid
[nl]
N=56
H1=75?m H1=100?m H1=125?m H1=150?m
10?1 100 101 102
4
6
8
10
12
14
16
18
20
?
V liquid
[nl]
N=30
H1=75?m H1=100?m H1=125?m H1=150?m
Figure 4.25: Variations of liquid charge with channel depth for different values of ?.
111
10?1 100 101 102
4
5
6
7
8
9
10
11
12
?
%V
MHP
N=60
H1=75?m
H1=100?m
H1=125?m
H1=150?m
10?1 100 101 102
4
5
6
7
8
9
10
11
12
?
%V
MHP
N=58
H1=75?m
H1=100?m
H1=125?m
H1=150?m
10?1 100 101 102
4
5
6
7
8
9
10
11
12
?
%V
MHP
N=56
H1=75?m
H1=100?m
H1=125?m
H1=150?m
10?1 100 101 102
4
5
6
7
8
9
10
11
12
?
%V
MHP
N=30
H1=75?m
H1=100?m
H1=125?m
H1=150?m
Figure 4.26: Percent of volume charge corresponding to maximum heat transport capacity
of the pipe.
112
10?1 100 101 102
4
6
8
10
12
14
16
?
q r[%]
N=60
H1=75?m H1=100?m H1=125?m H1=150?m
10?1 100 101 102
4
6
8
10
12
14
16
?
q r[%]
N=58
H1=75?m H1=100?m H1=125?m H1=150?m
10?1 100 101 102
4
6
8
10
12
14
16
?
q r[%]
N=56
H1=75?m H1=100?m H1=125?m H1=150?m
10?1 100 101 102
4
6
8
10
12
14
16
?
q r[%]
N=30
H1=75?m H1=100?m H1=125?m H1=150?m
Figure 4.27: Variations of MHP heat transport capacity with channel depth for different
values of ?.
113
Bibliography
[1] O.M. Alifanov. Inverse Heat Transfer Problems. Springer, 1994.
[2] K.K. Ambatipudi and M.M. Rahman. Analysis of conjugate heat transfer in microchan-
nel heat sinks. Numerical Heat Transfer, Part A, 37:711?731, 2000.
[3] K. Atkinson and W. Han. Theoretical Numerical Analysis- A Functional Analysis
Framework. Springer, 2001.
[4] M. Le Berre, G. Pandraud, P. Morfouli, and M. Lallemand. The performance of
micro heat pipes measured by integrated sensors. Journal of Micromechanics and
Microengineering, 16:1047?1050, 2006.
[5] R. Chein and G. Huang. Analysis of microchannel heat sink performance using nanoflu-
ids. Applied Thermal Engineering, 25(17-18):3104?3114, December 2005.
[6] E. K. P. Chong and S. Z?ak. An in introduction to optimization. Wiley, 2nd edition,
2001.
[7] M.J. Cola?co and H.R.B. Orlande. Inverse natural convection problem of simultaneous
estimation of two boundary heat fluxes in irregular cavities. International Journal of
Heat and Mass Transfer, 47(6-7):1201?1215, March 2004.
[8] T.P. Cotter. Principles and prospects of micro heat pipes. Proceedings 5th International
Heat Pipe Conference, Tsukuba, Japan:328?335, 1984.
[9] J.R. Culham, M.M. Yovanovich, and T.F. Lemczyk. Thermal characterization of elec-
tronic packages using three-dimensional Fourier series solutions. Journal of Electronic
Packaging, 122:233?239, September 2000.
[10] G. Dalaroz?ee. Introduction to reliability. Microelectronic Engineering, 49(1-2):3?10,
November 1999.
[11] M. Fabbri and V.K. Dhir. Optimized heat transfer for high power electronic cooling
using arrays of microjets. ASME Journal of Heat Transfer, 127:760?769, July 2005.
[12] A. Faghri. Heat Pipe Science and Technology. Taylor & Francis, 1995.
[13] C. Gillot, Y. Avenas, N. Cezac, G. Poupon, C. Schaeffer, and E. Fournier. Silicon heat
pipes used as thermal spreaders. IEEE Transactions on Components and Packaging
Technologies, 26(2):332?339, June 2003.
[14] C. Gillot, L. Meysenc, C. Schaeffer, and A. Bricard. Integrated single and two-phase
micro heat sinks under igbt chips. IEEE Transactions on Components and Packaging
Technologies, 22(3):384?389, September 1999.
114
[15] M. Groll, M. Schneider, V. Sartre, M.C. Zaghdoudi, and M. Lallemand. Thermal con-
trol of electronic equipment by heat pipes. Revue G?en?erale de Thermique, 37(5):323?
352, May 1998.
[16] M.D. Gunzburger and H.-C. Lee. Analysis and approximation of optimal control prob-
lems for first-order elliptic systems in three dimensions. Applied Mathematics and
Computation, 100(1):49?70, April 1999.
[17] J.M. Ha and G.P. Peterson. The interline heat transfer of evaporating thin films along
a micro grooved surface. ASME Journal of Heat Transfer, 118:747?755, August 1996.
[18] I. Hassan, P. Phutthavong, and Abdelgawad. Microchannel heat sinks: An overview
of the state-of-the-art. Microscale Thermophysical Engineering, 8:183?205, 2004.
[19] S. Hingorani, C.J. Fahrner, D.W. Mackowski, J.S. Goodling, and R.C. Jaeger. Optimal
sizing of planar thermal spreaders. ASME Journal of Heat Transfer, 116(2):296?301,
May 1994.
[20] C.-H. Huang and B.-H. Chao. An inverse geometry problem in identifying irregular
boundaryconfigurations. International Journal of Heat and Mass Transfer, 40(9):2045?
2053, June 1997.
[21] C.-H. Huang and H.-M. Chen. Inverse geometry problem of identifying growth of
boundary shapes in a multiple region domain. Numerical Heat Transfer, Part A,
35(4):435?450, March 1999.
[22] C.-H. Huang and J.-H. Hsiao. A non-linear fin design problem in determining the
optimum shape of spine and longitudinal fins. Communications in Numerical Methods
in Engineering, 19:111?124, 2003.
[23] C.-H. Huang and C.-C. Shih. A shape identification problem in estimating simulta-
neously two interfacial configurations in a multiple region domain. Applied Thermal
Engineering, 26(1):77?88, January 2006.
[24] C.-H. Huang and C.-Y. Yeh. An inverse problem in simultaneous estimating the Biot
numbers of heat and moisture transfer for a porous material. International Journal of
Heat and Mass Transfer, 45(23):4643?4653, November 2002.
[25] C.-H. Huang and C.-Y. Yeh. An optimal control algorithm for entrance concurrent flow
problems. International Journal of Heat and Mass Transfer, 46(6):1013?1027, March
2003.
[26] U. Imke. Porous media simplified simulation of single- and two- phase flow heat transfer
in micro-channel heat exchangers. Chemical Engineering Journal, 101(1-3):295?302,
August 2004.
115
[27] M. Ivanova, C. Schaeffer, Y. Avenas, A. Lai, and C. Gillot. Realization and thermal
analysis of silicon thermal spreaders used in power electronics cooling. IEEE Interna-
tional Conference on Industrial Technology, 2:1124?1129, 10-12 December 2003.
[28] D.S. Kercher, J.-B. Lee, O. Brand, M.G. Allen, and A. Glezer. Microjet cooling de-
vices for thermal management of electronics. IEEE Transactions on Components and
Packaging Technologies, 26(2):359?366, June 2003.
[29] D. Khrustalev and A. Faghri. Heat transfer during evaporation on capillary-grooved
structures of heat pipes. ASME Journal of Heat Transfer, 117:740?747, August 1995.
[30] D. Khrustalev and A. Faghri. Fluid flow effects in evaporation from liquid-vapor menis-
cus. ASME Journal of Heat Transfer, 118:725?730, 1996.
[31] D. Khrustalev and A. Faghri. Coupled liquid anf vapor flow in miniature passages with
micro grooves. ASME Journal of Heat Transfer, 121:729?733, 1999.
[32] R. W. Knight, D.J. Hall, J.S. Goodling, and R.C. Jaeger. Heat sink optimization
with application to microchannels. IEEE Transactions on Components, Hybrids, and
Manufacturing Technology, 15(5):832?842, October 1992.
[33] R.W. Knight, J.S. Goodling, and B.E. Gross. Optimal thermal design of air cooled
forced convection finned heat sinks- experimantal verification. IEEE Transactions on
Components, Hybrids, and Manufacturing Technology, 15(5):754?760, October 1992.
[34] A. Lai, C. Gillot, M. Ivanova, Y. Avenas, C. Louis, C. Schaeffer, and E. Fournier.
Thermal characterization of flat silicon heat pipes. 20th IEEE Semiconductor Thermal
Measurement and Management Symposium, pages 21?25, 9-11 March 2004.
[35] P. Lall, M. Pecht, and E. B. Hakim. Characterization of functional relationship between
temperature and microelectronic reliability. Microelectronics and Reliability, 35(3):377?
402, March 1995.
[36] S. Launay, V. Sartre, and M. Lallemand. Experimental study on silicon micro heat
pipe arrays. Applied Thermal Engineering, 24(2-3):233?243, February 2004.
[37] S. Lenhart and D.G. Wilson. Optimal control of heat transfer problem with convec-
tive boundary condition journal of optimization theory and applications. Journal of
Optimization Theory and Applications, 79(3):581?597, December 1993.
[38] H.-Y. Li and W.-M. Yan. Identification of wall heat flux for turbulent forced convection
by inverse analysis. International Journal of Heat and Mass Transfer, 46(6):1041?1048,
March 2003.
[39] J. Li, G. P. Peterson, and P. Cheng. Three-dimensional analysis of heat transfer in
a micro-heat sink with single phase flow. International Journal of Heat and Mass
Transfer, 47(19-20):4215?4231, September 2004.
116
[40] J.L. Lions. Optimal Control Systems Governed by Partial Differential Equations.
Springer, 1971.
[41] J.P. Longtin, B. Badran, and F. M. Gerner. A one-dimensional model of a micro heat
pipe during steady-state operation. ASME Journal of Heat Transfer, 116:709?715,
August 1995.
[42] J. C. De los Reyes. A primal-dual active set method for bilaterally control constrained
optimal control of the navier-stokes equations. Numerical Functional Analysis and
Optimization, 25(1):657?683, 2005.
[43] T. J. Martin and G. S. Dulikravich. Inverse determination of steady heat convection
coefficient distributions. ASME Journal of Heat Transfer, 120:328?334, May 1998.
[44] O.S. Nadgauda. Fabrication, filling, sealing and testing of micro heat pipes, Master
Thesis, Auburn University, December 2006.
[45] H.M. Park and W.J. Lee. A new numerical method for the boundary optimal con-
trol problems of the heat conduction equation. International Journal for Numerical
Methods in Engineering, 53(7):1593?1613, March 2002.
[46] H.M. Park and H.J. Shin. Empirical reduction of modes for the shape identification
problems of heat conduction systems. Computer Methods in Applied Mechanics and
Engineering, 192(15):1893?1908, April 2003.
[47] K. Park and K.-S. Lee. Flow and heat transfer characteristics of the evaporating
extended meniscus in a micro-capillary channel. International Journal of Heat and
Mass Transfer, 46(24):4587?4594, November 2003.
[48] Y. Peles, A. Kosar, C. Mishra, C.-J. Kuo, and B. Schneider. Forced convective heat
transfer across a pin fin micro heat sink. International Journal of Heat and Mass
Transfer, 48(17):3615?3627, August 2005.
[49] X. F. Peng and G. P. Peterson. The effect of thermofluid and geometrical parameters
on convection of liquids through rectangular microchannels. International Journal of
Heat and Mass Transfer, 38(4):755?758, March 1995.
[50] X. F. Peng and G. P. Peterson. Convective heat transfer and flow friction for water
flow in microchannel structures. International Journal of Heat and Mass Transfer,
39(12):2599?2608, August 1996.
[51] X. F. Peng, B. X. Wang, G. P. Peterson, and H. B. Ma. Experimental investigation
of heat transfer in flat plates with rectangular microchannels. International Journal of
Heat and Mass Transfer, 38(1):127?137, January 1995.
117
[52] P. Razelos and R.N. Krikkis. On the optimum thermal design of individual longitudinal
finswith rectangular profile. International Communications in Heat and Mass Transfer,
30(3):349?358, 2003.
[53] J. H. Ryu, D. H. Choi, and S. J. Kim. Numerical optimization of the thermal perfor-
mance of a microchannel heat sink. International Journal of Heat and Mass Transfer,
45(13):2823?2827, June 2002.
[54] J. H. Ryu, D. H. Choi, and S. J. Kim. Three-dimensional numerical optimization of
a manifold microchannel heat sink. International Journal of Heat and Mass Transfer,
46(9):1553?1562, April 2003.
[55] A.K. Saha and S. Acharya. Parametric study of unsteady flow and heat transfer in a
pin-fin heat exchanger. International Journal of Heat and Mass Transfer, 46(20):3815?
3830, September 2003.
[56] O. Salmela. The effect of introducing increased-reliability-risk electronic components
into 3rd generation telecommunications systems. Reliability Engineering & System
Safety, 89(2):208?218, August 2005.
[57] V. Sartre, M.C. Zaghdoudi, and M. Lallemand. Effect of interfacial phenomena on
evaporative heat transfer in micro heat pipes. International Journal of Thermal Sci-
ences, 39:498?504, 2000.
[58] M. Silieti, E. Divo, and A.J. Kassab. An inverse boundary element method/genetic
algorithm based approach for retrieval of multi-dimensional heat transfer coefficients
within film cooling holes/slots. Inverse Problems in Science and Engineering, 13(1):79?
98, February 2005.
[59] F. Simionescu and D.K. Harris. Estimation of the heat transfer coefficient of
a microchannel heat sink using an optimal control technique. In Proceedings of
ASME ICNMM2006 4th International Conference on Nanochannels, Microchannels
and Minichannels, Limerick, Ireland, June 19-21, 2006.
[60] F. Simionescu, A.J. Meir, and D.K. Harris. Approximation of an optimal convective
heat transfer coefficient. Optimal Control Applications and Methods, 27:237?253, 2006.
[61] P.C. Stephan and C.A. Busse. Analysis of the heat transfer coefficient of grooved heat
pipe evaporator walls. International Journal of Heat and Mass Transfer, 35(2):383?391,
1992.
[62] B. Suman, S. De, and S. DasGupta. A model of the capillary limit of a micro heat pipe
and prediction of the dry-out length. International Journal of Heat and Fluid Flow,
26:495?505, 2005.
118
[63] B. Suman and N. Hoda. Effect of variations in thermophysical properties and design
parameters on the performance of a V-shaped micro grooved heat pipe. International
Journal of Heat and Mass Transfer, 48:2090?2101, 2005.
[64] L.W. Swanson and G.P. Peterson. The interfacial thermodynamics of micro heat pipes.
ASME Journal of Heat Transfer, 117:195?201, February 1995.
[65] S.K. Thomas, R.C. Lykins, and Kirk L. Yerkes. Fully developed laminar flow in trape-
zoidal grooves with shear stress at the liquidvapor interface. International Journal of
Heat and Mass Transfer, 44(18):3397?3412, September 2001.
[66] C.L. Tien, A. Majumdar, and F. Gerner. Microscale Energy Transport. Taylor &
Francis, 1998.
[67] K.C. Toh, X.Y. Chen, and J.C. Chai. Numerical computation of fluid flow and
heat transfer in microchannels. International Journal of Heat and Mass Transfer,
45(26):5133?5141, December 2002.
[68] D. B. Tuckerman and F. W. Pease. High performance heat sinking for VLSI. IEEE
Electron Device Letters, 2(5):126?129, May 1981.
[69] K. Vafai and L. Zhu. Analysis of two-layered micro-channel heat sink concept in
electronic cooling. International Journal of Heat and Mass Transfer, 42(12):2287?2397,
June 1999.
[70] C.-Y. Wang, M. Groll, S. R?osler, and C.-J. Tu. Porous medium model for two-phase
flow in mini channels with applications to micro heat pipes. Heat Recovery Systems &
CHP, 14(4):377?389, 1994.
[71] J.J. Wei and H. Honda. Effects of fin geometry on boiling heat transfer from silicon
chips with micro-pin-fins immersed in FC-72. International Journal of Heat and Mass
Transfer, 46(21):4059?4070, October 2003.
[72] A. Weisberg, H. H. Bau, and J. N. Zemel. Analysis of micro-channels for integrated
cooling. International Journal of Heat and Mass Transfer, 35(10):2465?2474, 1992.
[73] H. Y. Wu and P. Cheng. An experimental study of convective heat transfer in silicon
microchannels with different surface conditions. International Journal of Heat and
Mass Transfer, 46(14):2547?2556, July 2003.
[74] R.-H. Yeh and M. Chang. Optimum longitudinal convective fin arrays. International
Communications in Heat and Mass Transfer, 22(3):445?460, 1995.
[75] E. P. Zafiropoulos and E. N. Dialynas. Reliability and cost optimization of electronic
devices considering the component failure rate uncertainty. Reliability Engineering &
System Safety, 84(3):271?284, June 2004.
119
[76] C. Y. Zhao and T. J. Lu. Analysis of microchannel heat sinks for electronics cooling.
International Journal of Heat and Mass Transfer, 45(24):4857?4869, November 2002.
120