Complete analytical model of a loop heat pipe with a

flat evaporator
B. Siedel, V. Sartre, F. Lefèvre

To cite this version:

B. Siedel, V. Sartre, F. Lefèvre. Complete analytical model of a loop heat pipe with a flat evaporator.
International Journal of Thermal Sciences, 2015, 89, pp.372-386. �10.1016/j.ijthermalsci.2014.11.014�.
�hal-01286580�

HAL Id: hal-01286580

https://hal.science/hal-01286580
Submitted on 11 Mar 2016

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Complete analytical model of a loop heat pipe with a flat evaporator

Benjamin Siedel, Valérie Sartre∗, Frédéric Lefèvre

Université de Lyon, CNRS
INSA-Lyon, CETHIL UMR5008, F-69621, Villeurbanne, France
Université Lyon 1, F-69622, Villeurbanne, France

Abstract
A steady-state analytical model has been developed to determine the thermohydraulic behaviour of a loop heat pipe
with a flat evaporator. Its main originality lies in the combination of energy balance equations for each component of the
system with 2D analytical solutions for the temperature field in the evaporator. Based on Fourier series expansion, heat
transfer in the wick as well as in the evaporator casing are accurately modelled, enabling a thorough consideration of the
parasitic heat fluxes. The model is based on the thermal contact resistance between the wick and the casing, the thermal
conductivity of the wick and the accommodation coefficient. This analytical method offers a simple solution that can
be implemented in LHP design analysis without the need of large computational resources. A sensitivity analysis has
been carried out to evaluate the influence of several parameters on the LHP behaviour. The results show that the main
parameters of the model are independent. Therefore, they could be experimentally determined using an appropriate test
bench with only few temperature measurements. The model has been validated with a set of experimental data from
the literature. A good agreement is found between the theoretical and the experimental results.
Keywords: Loop Heat Pipe, Analytical, Modeling, Sensitivity analysis

1. Introduction the advantage to provide a solution without an excessive

computing time and that can be easily implemented. Some
Loop Heat Pipes (LHP) are two-phase cooling systems analytical models of LHP are found in research works.
able to passively transport high amounts of heat over dis- Maydanik et al.[7] (cited by Launay et al.[8]) suggested
tances up to several meters. Developed in the 1970’s, these a closed form solution of a LHP analytical model estab-
devices have proven their reliability in many spatial appli- lishing the energy balance in the reservoir and the pres-
cations and are now candidates for terrestrial cooling so- sure balance in the whole system. Furukawa[9] developed
lutions. Indeed, their specific design offers robustness and a complete analytical model of a LHP able to enhance the
flexibility for a wide variety of practical applications [1, 2]. sizing of the system and to study the influence of many
As a consequence, many efforts have been dedicated geometrical parameters on the loop operation. However,
to understand their operation in order to optimize their the model requires the operating temperature as input pa-
design. A LHP is a complex thermal system including rameter. Yet in most cases, the evaporator temperature
an evaporator connected to the heat source, a condenser is the main expected output of a LHP model. Launay
to dissipate the heat load and vapour and liquid lines et al.[8] proposed closed-form expressions of the operat-
to transport the working fluid between both components. ing temperature of the LHP, for both the variable and the
Coupled thermo-hydraulic phenomena govern the LHP be- fixed conductance mode. Their model is based on energy
haviour and need to be understood to enable a correct sys- balance equations on each system component and on ther-
tem designing. Parameters such as evaporation efficiency, modynamic equations. The thermal links in the reservoir
heat losses to the ambient and parasitic heat flux in the are defined as equivalent thermal resistances. Their solu-
evaporator as well as condensation heat transfer can be of tion is a useful tool in the LHP design. However, the heat
great influence on the loop operation. transfer in the evaporator is not accurately determined and
A lot of papers concerning LHP complete modelling can has to be adjusted in accordance with experimental data.
be found in the literature [3–6]. However, these numerical The purpose of the present study is to present a com-
models often imply complex algorithms and large compu- plete analytical model of a Loop Heat Pipe accurately tak-
tational resources which are not always available for up- ing into account heat and mass transfer in the evaporator
stream pre-design applications. An analytical model offers structure. The model is developed for a flat disk-shaped
evaporator geometry. This model is based on the analyti-
∗ Corresponding
cal study of Launay et al.[8]. However, the present study
author
Email address: valerie.sartre@insa-lyon.fr (Valérie Sartre) rests upon two 2-D analytical solutions describing the tem-

Preprint submitted to International Journal of Thermal Sciences January 10, 2014

perature field in the wick and the evaporator casing. Both Condenser
solutions enable to determine parasitic heat losses through Evaporator/Reservoir two-phase
the wick and the evaporator body, the sensible heat given
to the liquid flowing through the porous structure as well Q
as the heat dissipated by evaporation at the wick-groove Body in Wick
interface. A similar approach was implemented in the case Text Qb Qw Qev Tv Tsink
of conventional heat pipes by Lefèvre and Lallemand [10] Qext,e Te Twe
and later extended by Lips and Lefèvre [11]. These fea-
tures, coupled with energy balances and thermodynamic
relationships in the rest of the LHP, give a simple solution
for the operating temperature. An iterative procedure is
implemented to calculate the two-phase length in the con-
denser.
Text
∆P
2. Model description Qext,r Tr Qsub

2.1. Analytical Model of the LHP

The thermal state of the complete LHP can be deter-
mined using energy balance equations and thermodynamic
relationships. Figure 1 presents the operating principle of Tr,in
the LHP and the links between its components.

Liquid line
The total heat load to be dissipated by the evaporator Thermodynamic
∆P relationship
Qin is conducted through the wick or through the evapo- Text Fluid thermal
rator body so that:
resistance
Qin = Qw + Qb (1) Conduction
The wick is assumed to be fully saturated with liquid. The Tc,o Convection
thermal heat flux Qw is transversally conducted through
subcooling
Condenser

the evaporator wall, the wall-wick interface and then di- Phase-change
vides up: a part Qev is evaporated at the wick-groove inter- Tsink
face whereas the rest is dissipated by conduction and con- Contact
vection with the liquid flowing through the porous struc-
ture and with the liquid in the reservoir. Qb is conducted
longitudinally through the evaporator wall to the reservoir Tv
and a part of it, Qext,e , is given by convection to the am-
bient. Both the heat flux through the wick Qw and the
heat flux conducted through the evaporator casing Qb are Figure 1: LHP schematic nodal network
functions of the reservoir, the groove, the wick and the
evaporator temperatures Tr , Tv , Twe and Te . The same
dependence applies for Qev and Qext,e :
reservoir whereas the rest is dissipated by convection due
Qw = f (Tr , Tv , Te ) (2) to the liquid flow inside the porous structure. At the inter-
Qev = ṁl hlv = f (Tr , Tv , Twe ) (3) face between the wick and the evaporator envelope, there
is a temperature gap Te -Twe due to a contact resistance Rc
Qb = f (Tr , Te ) (4)
defined as:
Qext,e = f (Tr , Te ) (5) Te − Twe
Rc = Sc (7)
Qw
Qext,e is also a function of Text , which is a given data of where Sc is the contact surface between the wick and the
our model. As a result, the heat load Qin can also be evaporator body and Te and Twe are the temperatures on
expressed as a function of these four temperatures: the envelope side and on the wick side, respectively.
A global heat balance on the evaporator/reservoir gives
Qin = f (Tr , Tv , Twe , Te ) (6)
the following equation:
An analytical expression of Qin will be derived in subsec- Qin = Qev + Qsen + Qsub + Qext,e + Qext,r (8)
tions 2.2 and 2.3. The part of Qw that is not dissipated by
evaporation is the transversal parasitic heat flux. Part of where Qsen is the sensible heat given to the liquid, Qsub
this flux is conducted through the wick and released to the is the subcooling due to the liquid entering the reservoir
2
and Qext,r is the heat flux dissipated to the ambient by the are calculated considering a convective heat transfer hsink
reservoir. The determination of Qsen and Qsub leads to: and hext , respectively:

Qsen = ṁl cp,l (Tv − Tr ) (9) Tc,o = Tsink + (Tv − Tsink ) (16)
 
Qsub = ṁl cp,l (Tr − Tr,in ) (10) −πDc,i (Lc − L2φ )
× exp
ṁl cp,l (1/hl + Dc,i /(hsink Dc,o ))
where Tr,in is the temperature of the liquid coming from
the condenser and flowing back to the reservoir. where Tc,o is the temperature of the liquid at the condenser
To evaluate the heat transfer given by the reservoir to the outlet, Lc and L2φ are the lengths of the condenser and
ambient, it is assumed that its surface is at a uniform of the two-phase region respectively. Dc,i and Dc,o are the
temperature equal to Tr . Heat transfer with the ambient inner and outer condenser diameters. In a similar way, the
Qext,r is then approximated by: liquid line heat balance is:
Tr,in = Text + (Tc,o − Text ) (17)
Qext,r = hext Sr (Tr − Text ) (11)  
−πDl,i Ll
× exp
where Sr is the total external surface of the reservoir. A ṁl cp,l (1/hl + Dl,i /(hext Dl,o ))
single expression of Qin can be derived by combining equa-
where Dl,i and Dl,o are the inner and outer diameters of
tions (3), (5) and (7) to (11):
the liquid line and Ll its length. Additionally, in the con-
Qin = f (Tr , Tv , Twe , Te , ṁl , Tr,in ) (12) denser, the heat exchange with the heat sink in the two-
phase zone is equal to the latent heat of the condensing
In some cases, the reservoir can be full of liquid. This phe- vapour:
nomenon has been extensively studied by Adoni et al.[12]. 1
In the present model, the existence of a two-phase equilib- ṁl hlv = 1 1 πL2φ (Tv − Tsink ) (18)
+
rium in the reservoir is assumed. As shown by Launay et hcond Dc,i hsink Dc,o

al. [8], a thermodynamic relationship exists between the Equations (3), (7), (8), (12), (15) to (17) and (18) form
saturation temperature inside the grooves and the one at a set of 8 independent equations with 8 unknowns: Twe ,
the liquid-vapour interface in the reservoir: Te , Tr , Tr,in , Tc,o , Tv , ṁl and L2φ . Their solution leads
  to the determination of the complete thermal state of the
∂T
∆T = Tv − Tr = (∆Pv + ∆Pl − ρl g∆H) (13) LHP. The values of Qb , Qw and Qev (equations 2-4) can
∂P
be calculated using several methods. In the large majority
where ρl is the liquid density and ∆H is the elevation of of LHP models from the literature, heat transfer in the
the condenser compared to the evaporator. The slope of evaporator is simply described using equivalent thermal
the pressure-temperature saturation curve is given by the resistances based on the geometrical characteristics and
Clausius-Clapeyron equation: the thermophysical properties of the evaporator. However,
such a method does not take into account adequately the
∂T T (1/ρv − 1/ρl ) heat flux in the wick and the evaporator body since the
= (14)
∂P hlv determination of the thermal resistances requires a 2D or
3D approach. Thus, in the present study, an accurate
The model can also cope with non-condensable gases
thermal analysis of the evaporator has been conducted,
(NCG), which can be generated for various reasons in the
based on an analytical approach.
LHP and accumulate in the reservoir, as explained by
Singh et al. [13]. We assume that in operating condi- 2.2. Analytical Thermal Model of the Wick
tions, these NCG are drained to the reservoir. In order to
The first part of the analytical model of the evaporator
take into account the overpressure generated by the NCG
deals with heat and mass transfer inside the porous struc-
PNCG , equation (13) is modified:
ture. As shown in Figure 2, a part of the porous wick is

∂T
 modelled, bordered on one side by the liquid bulk of the
∆T = Tv − Tr = (∆Pv + ∆Pl − ρl g∆H + PNCG ) reservoir and by half of a fin and half of a groove on the
∂P
(15) other side.
The vapour line is considered adiabatic so that the The 2-D stationary heat equation in the wick is expressed
vapour enters the condenser with a temperature equal to as:
∂ 2 Tw ∂ 2 Tw
Tv . Furthermore, the condensation temperature Tc and + =0 (19)
the vapour temperature Tv are linked with the thermody- ∂x2 ∂y 2
namic relationship (14). Since the pressure drops in the where Tw is the temperature of the porous structure and x
vapour line are low, it is assumed that condensation occurs and y are the axis coordinates (Fig. 3). A non-dimensional
at temperature Tv (Tc ≈ Tv ). temperature is defined as:
In the part of the condenser where liquid subcooling oc- λeff (Tw − Tr )
curs, heat transfer with the heat sink and with the ambient Tw∗ = (20)
φ0 b
3
 Thus equation (23) becomes:
Reservoir vapour
∞
X
liquid Tw∗ (X, Y ) = Am (Y )cos(mπX) (25)
m=0
– – – b– – – – – – – – – – – – –
– – – – – – – – – – – – – – – – Combining equations (21) and (25) yields:
– – – – – – – – – – – – – – – –
– – – – – – – – – – – – – – – – ∂ 2 Am (Y )
– wick − (mπB)2 Am (Y ) +
– – – – – – – – – – – – – – – –
∂Y 2
=0 (26)
– – – – – – – – – – – – – – – –
– – – – – – – – – – – – – – – – Solving the previous differential equation leads to:
– – – – – – – – –
a – – – – – –
groove A0 (Y ) = A01 Y + A02 if m = 0 (27)
mπBY
wall Am (Y ) = Am1 e + Am2 e −mπBY
otherwise (28)

The boundary condition for Y = 0 is a set temperature

Figure 2: Evaporator cross-section profile corresponding to the evaporator temperature Twe
for one side (x 6 a0 ) and to the groove temperature Tv for
the other side (x > a1 ). At the junction between the wall
and the groove, the temperature singularity is treated by
and leads to: considering a linear variation of the temperature between
∂ 2 Tw∗ 1 ∂ 2 Tw∗ two fictive points a0 and a1 . The slope of the temperature
+ =0 (21)
∂X 2 B 2 ∂Y 2 gradient and thus the position of these points depend on
the maximum evaporation rate at the liquid-vapour inter-
with
face that can be calculated using the kinetic gas theory,
x y b
X= , Y = and B= (22) which gives a heat transfer coefficient hev [14]:
a b a
−0.5 
2aev ρv h2lv
 
where a and b are the lengths of the modelled region in the 2πRTsat Psat
hev = 1− (29)
x and y directions respectively, λeff is the wick effective 2 − aev Tsat M 2ρv hlv
thermal conductivity and φ0 is an arbitrary heat flux.
where aev is the accommodation coefficient. Such an ap-
proach was already used in [15] to cope numerically with
y the temperature singularity at the triple line. In the case
of the evaporation of a thin liquid film, the accommodation
Tr coefficient is defined as the ratio of the actual evaporation
br
r r rate to a theoretical maximal phase change rate. A coef-
r r
r
r
r r ficient equal to unity describes perfect evaporation while
r r
r
r
r r a lower value represents incomplete evaporation. In the
r r
r
r
r r case of water, values varying from 0.01 to 1 are suggested
r r
r
r
r r in the literature [16]. We assume the following relationship
r r
r
r
r r between hev and the position of a0 and a1 :
r r
r
r r
r r λeff
0 Twe a0 a1 Tv a x hev = (30)
a1 − a0
Figure 3: Schematic of the wick model The boundary condition for Y = 0 is then defined as:
λeff (Twe − Tr ) a0


 if 0 < X 6
φ0 b a

A general expression of the non-dimensional temperature




field expanded in 2-D Fourier series is given by:
  
λ a − aX


 eff (Tv − Tr ) + (Twe − Tv ) 1


∞
X Tw∗ (X, 0) = φ0 b a0
a1 − a0
a1
Tw∗ (X, Y ) = Am (Y )cos(mπX) + Bm (Y )sin(mπX)

 if < X <
a a


m=0



(23)

λ (T − Tr ) a1

 eff v

if 6X<1

For X = 0 and X = 1, an adiabatic boundary condition is

φ0 b a
considered based on a symmetry hypothesis: (31)
For Y = 1, the reservoir temperature is set:
∂Tw∗ ∂Tw∗
= =0 (24)
∂X X=0 ∂X X=1
Tw∗ (X, 1) = 0 (32)
4
The non-dimensional temperature field is then: The previously described model enables to calculate the
heat flux through the wick Qw and the heat dissipated by
Tw∗ (X, Y ) = A01 Y + A02 evaporation Qev :
(33)
∞
X
mπBY


+ Am1 e + Am2 e −mπBY
cos(mπX) a0 +a1
Sw ∂Tt
Z 2
m=1 Qw = −λeffdx (44)
a 0 ∂y y=0
with   
a0 + a1 3qb
  = −Sw φ0 A01 −
λeff a0 + a1 2a 2φ0
A01 = − Tv − Tr + (Twe − Tv ) (34)
φ0 b 2a ∞  #
X a0 + a1
+ B (Am1 − Am2 ) sin mπ
 
λeff a0 + a1 2a
A02 = Tv − Tr + (Twe − Tv ) (35) m=1
φ0 b 2a
cos mπ aa0 − cos mπ aa1



λeff Twe − Tv a
Am1 =2
φ0 b m2 π 2 a1 − a0 1 − e2mπB Sw a
∂Tt
Z
(36) Qev = λeff dx (45)
a a0 +a1
2
∂y y=0
a0 a1



λeff Twe − Tv a cos mπ a − cos mπ a 
a0 + a1

3qb

Am2 =2 = Sw φ0 1 − A01 −
φ0 b m2 π 2 a1 − a0 1 − e−2mπB 2a 2φ0
(37) ∞  #
X a0 + a1
− B (Am1 − Am2 ) sin mπ
In equation (33), the liquid flow inside the porous 2a
m=1
structure is not taken into account. This flow is two-
dimensional. However, in the present study, we assume
a 1D flow inside the wick, considering a homogeneous vol- 2.3. Analytical Thermal Model of the Evaporator Body
umetric source inside the wick q defined as:
The model of the wick is sufficient if the heat losses
4ṁl cp,l (Tv − Tr ) through the evaporator body are negligible. Neverthe-
q=− 2b
(38)
πDw less, the body is usually made of a high conductive ma-
where ṁl is the total liquid mass flow rate in the wick terial. Therefore, it is generally necessary to estimate the
and cp,l is the specific heat of the liquid. This assumption heat transferred by conduction from the evaporator to the
respects the energy balance and enables to derive an ana- reservoir and to the ambient. A second analytical model
lytical expression for the influence of the liquid flow inside is developed to describe the evaporator casing. As shown
the wick. The superposition principle enables to add a in Figure 4, the evaporator wall is “unwrapped” and mod-
simple model of heat transfer in the wick with a homoge- elled as a rectangular domain. At x = 0 and x = c, an
neous source to the previous developed analysis. It has to adiabatic condition is assumed due to the symmetry. At
be noted that a 2D approach of heat and mass transfer in y = 0, a convective heat transfer with the ambient is taken
a porous wick was presented by Cao and Faghri [17], but into account on the whole external surface, including the
the model is not entirely analytical. The heat equation in heating section (0 6 x 6 c0 ), where a heat flux φin is also
the wick is thus: applied. The inner part of the body is highly influenced by
∂ 2 Tc q the reservoir temperature. The rectangular shape chosen
= (39) to represents the evaporator body does not consider the
∂y 2 λeff
effect of the grooves on the distortion of the temperature
and the boundary conditions are:
field between the base plate and the wick. We assume,
q therefore, a set temperature profile, with the evaporator
Tc (y = 0) − Tc (y = b) = b2 (40)
λeff temperature Te on one side (0 6 x 6 c0 ) and the reservoir
Tc (y = b) = Tr (41) temperature Tr on the other side (c1 6 x 6 c) with a linear
profile in-between.
This leads to the non-dimensional temperature field im- The sum of two solutions is necessary to take into ac-
plied by the liquid flow inside the wick: count the complete set of boundary conditions. The first
one corresponds to a set temperature profile in the inner
qb
Y 2 − 3Y + 2


Tc∗ = (42) part of the wall, in contact with the wick and the liquid-
2φ0 vapour bulk in the reservoir. Convective heat losses to the
This solution is added to the previous one to give a general ambient are assumed for the whole external surface of the
expression of the non-dimensional temperature field in the evaporator/reservoir. This hypothesis considers that the
wick, with consideration of 1-D convection: electronic component to cool down dissipates heat with the
LHP on one side and with the ambient on the other. The
Tt∗ = Tw∗ + Tc∗ (43) second solution adds the heat input to the evaporator.
5
 y
vapour dr Te Tr
r r
r
r
r r
r Tb1 (x, y) First solution r
r
r
r r
r
liquid Tr 0 c0 c c x
hext1
Text
c1
– – – – – – – –
– – – – – – – – –
– – – – –
wick – – –
– – – – – – – – –
=⇒ +
– – – – – – – – y

Qb
– – – – – – – – – Text
dr T =0
Te r r
r
hext r
r r
r Tb2 (x, y) Second solution r
r
r
r r
r
c0 0 c0 φ=0 c x
φin φin

Figure 4: Evaporator casing modelling

Using the same mathematical procedure as previously, The combination of equations (46) to (53) and the addition
the heat equation becomes: of both solutions lead to the non-dimensional temperature
field in the evaporator body:
∂ 2 Tb∗ 1 ∂ 2 Tb∗
+ =0 (46) Tb∗ (X, Y ) = C01 Y + C02
∂X 2 D2 ∂Y 2 (54)
∞
where:
X
mπDY


+ Cm1 e + Cm2 e −mπDY
cos(mπX)
m=1
x y d
λb (Tb − Tr )
X= ; Y = ; D= ; Tb∗ =
c d c φ0 d with
(47)  
The boundary conditions for the first solution Tb1
∗
are: Bib λb c0 + c1
C01 = (Te − Tr ) + (Tr − Text )
Bib + 1 φ0 d 2c
∂Tb1
∗
∂Tb1
∗
c0 φin
= =0 (symmetry) (48) − (55)
∂X X=0 ∂X X=1 c φ0
 
1 λb c0 + c1
∂Tb1
∗
hext C02 = (Te − Tr ) − Bib (Tr − Text )
= Bib Tb1
∗
+ (Tr − Text ) Bib + 1 φ0 d 2c
∂Y Y =0 φ0 c0 φin
(49) + (56)
dhext c φ0
with Bib =
λb cos mπ cc0 − cos mπ cc1



λb Te − Tr c
Cm1 = 2 2 2 mπD+Bib −mπD

λ (T − Tr ) c0 φ0 d m π c1 − c0 emπD + mπD−Bi e
 b e
 if 0 < X 6
b

 φ0 d c φin 1  c0  1

 −2 2 2
sin mπ (57)
φ0 m π D c 1 + e2mπD

Tb1
∗
(X, 1) = λb c1 − cX c0 c1
(Te − Tr ) if 6X<
cos mπ cc0 − cos mπ cc1



 φ0 d c1 − c0 c c λb Te − Tr c
Cm2 = 2


c1 φ0 d m2 π 2 c1 − c0 e−mπD + mπD−Bi
 b mπD
mπD+Bib e

0

if <X<1
c φin 1  c0  1
(50) +2 sin mπ (58)
2
φ0 m π D 2 c 1 + e−2mπD
whereas for the second solution Tb2
∗
:
The heat dissipated to the ambient in the evaporator
∂Tb2
∗
∂Tb2
∗ section of the body is:
= =0 (51)
∂X X=0 ∂X X=1 Z c0
 φ Qext,e = hext πx (Tb (x, 0) − Text ) dx (59)
in c0
∂Tb2∗ −
 if 0 < X 6 Z 0c1
= φ 0 c (52)
∂Y Y =0  c 0
+ hext πc0 (Tb (x, 0) − Text ) dx
0 if X > c0
c  c0 

φ0 d

Tb2
∗
(X, 1) = 0 (53) = hext πc0 c1 − Tr − Text + C02
2 λb
6
 ∞
φ0 d X c Parameter initialization, Qin (1)
+ hext π (Cm1 + Cm2 )
λb m=1 mπ
 c   mπc  c0  c1 
0
× −1 + cos + sin mπ Setting of the two-phase
mπ 2 mπ c length in the condenser
L2φ according to Eq. (18)
The heat transferred by thermal conduction through the
evaporator body to the reservoir section is calculated by Calculation of the thermophysical
integrating the Fourier’s law at y = c1 (Fig. 4). Thus, the properties and the pressure drops
total heat losses through the body is:
Calculation of the K-
d coefficients (See Appendix A)
∂Tb
Z
Qb = Qext,e + −πc0 λb dy (60)
0 ∂x x=c1
Determination of the vapour
∞ temperature Tv (Eq. B.1)
X  c1 
= Qext,e + dc0 φ0 π sin mπ
m=1
c
× Cm1 emπD − 1 − Cm2 e−mπD − 1




Is the energy balance in no
the condenser (Eq. 18)
satisfied?
2.4. Solving procedure
The solving procedure is presented in Figure 5. The set yes
of equations is not linear. Thus an iterative procedure is
Qin (i + 1) = Qin (i) + ∆Qin
used to solve it. After initialisation of the parameters of
the model, the two-phase length in the condenser L2φ is
set, according to the energy balance for a given heat input
Qin (18). L2φ has a major influence on the determination no
Qin (i + 1) > Qmax
of the temperature of the liquid entering in the reservoir
Tr,in (equations 16-17). Then, the thermophysical prop- yes
erties are calculated, as well as the pressure drops in the
transport lines. Plot of the operating curve
K-coefficients can be defined to reformat the expression
Figure 5: Solving algorithm flowchart
of Qw , Qev , Qb and Tr,in as functions of Tr , Tv , Te and
Twe :

Qw = K1 Tr + K2 Tv + K3 Twe using (34,36,37,44)

(61) the vapour and the liquid lines, depend on the fluid flow
regime. They are calculated as follows:
Qev = K4 Tr + K5 Tv + K6 Twe using (34,36,37,45)
(62)  2
Qext,e = K10 Tr + K11 Te + K12 using (56,57,58,59) f ṁ
∆P = L (66)
(63) 2ρD A
Qb = K7 Tr + K8 Te + K9 using (57,58,60,63)
(64) where A is the cross-section area of the tube. For a smooth
tube wall of diameter D and length L, the friction factor
Tr,in = K13 Tv + K14 using (16,17) (65) f is expressed by:
The detailed expression of these coefficients is presented 
in Appendix A. Equations (1,7-11, 15-17,61-65) are solved 64/Re
 if Re ≤ 2000
and give a second-order expression which enables to cal- f = 0.032 if 2000 < Re < 9150 (67)
culate the vapour temperature Tv (see Appendix B). 

0.316Re−0.25 if Re ≥ 9150
This procedure is iterated until the energy balance is
satisfied in the condenser (equation 18). The same method
is computed for each heat input Qin (i). Since the heat Equation (15) also includes the pressure of non-
transfer coefficient with the heat sink is generally much condensable gases PNCG . In order to take into account
lower than the condensation heat transfer coefficient, the the NCG, it is necessary to calculate the liquid level in the
condensation thermal resistance is neglected to simplify reservoir that depends on the heat load. We assume that
equation (18). the void fraction of the two-phase flow in the condenser is
Equation (15) is a function of the pressure drops ∆Pv 0.5 and that the vapour density is negligible compared to
and ∆Pl . These parameters, due to the friction forces in the liquid density. The liquid level el in the reservoir is
7
therefore expressed as: Extensive studies have been undertaken to develop mod-
! els able to predict accurately the effective thermal conduc-
1 mf 2
Ll πDl,i (Lc − 21 L2φ )πDc,i
2
tivity of a porous structure. Thus, many different corre-
el = − εSw ew − −
Sw ρl 4 4 lations can be found in the literature to determine the
(68) effective thermal conductivity of a porous material [21–24]
where mf is the total fluid charge in the system and ew is Table 2 presents the results obtained for a nickel wick of
the wick thickness. The total volume of NCG and vapour 75 % porosity saturated with water, using various correla-
in the reservoir Vv is equal to: tions. This thermophysical property depends not only on
the conductivity of the materials constituting the wick but
Vv = Sw (er − el ) (69) also on geometrical parameters such as the porosity, the
mean pore diameter and the pore size distribution. The
with er the thickness of the reservoir. Considering the non- results are very different according to the chosen correla-
condensable gases as ideal gases, their partial pressure is tion. This shows that this parameter is difficult to evalu-
calculated by: ate accurately. For the standard case, an effective thermal
mNCG RTr conductivity equal to 5 W · m=1 · K=1 is chosen, which cor-
PNCG = (70)
MNCG Vv responds to the wick properties of Singh et al.[19].

λeff (W · m=1 · K=1 )

where R is the ideal gas constant and mNCG and MNCG
are the NCG total mass and the molar mass, respectively. Correlation
To solve equations (16-18), the heat sink heat transfer Alexander 5.82
coefficient hsink is determined in accordance with the con- Chaudhary-Bandhari 4.06
denser design, whereas the heat transfer coefficient of the Maxwell 16.96
liquid hl is calculated assuming a laminar fully-developed Parallel scheme 23
flow and a constant Nusselt number N uD = 4.36 [18]. Zehner-Schlunder 1.75
Equations (11) and (17) depend on the heat transfer
coefficient with the ambient hext . It is given by the corre- Table 2: Effective thermal conductivity calculation with a water-
lation of Churchill and Chu [18] for free convection on the saturated nickel wick, 75% porosity
surface of an isothermal cylinder:
 2 Information from the literature concerning the thermal
1
λair  0.387RaD 6 contact between a porous structure and a solid base plate
hext = 0.60 +  (71) is very scarce. This parameter depends on many geomet-

8 
D 9
 27
1 + (0.559/P r) 16 rical and manufacturing characteristics. Choi et al.[20] in-
vestigated new techniques to enhance the thermal contact
This correlation is valid for Rayleigh numbers RaD lower conductance of evaporators in LHPs. Several channel de-
than 1012 . The determination of the effective thermal con- signs have been tested and compared. An estimation of the
ductivity of the porous structure will be discussed in the contact thermal resistance (including in practice thermal
next section. conduction in the base plate, in the wick and an evapora-
tion resistance) is given and values ranging from 7 · 10=5
to 3 · 10=4 m2 · K · W=1 have been obtained. Based on this
3. Results and discussion work, a constant contact resistance of 10=4 m2 · K · W=1 is
chosen for the standard case.
This section presents a sensitivity analysis to show the
Although the accommodation coefficient is a key pa-
influence of different parameters on the LHP performance.
rameter in the determination of the LHP temperature, its
The LHP geometry considered for this analysis is similar
value is very difficult to predict. In the literature, several
to standard systems used for electronic cooling applica-
studies are dedicated to the determination of this param-
tions. Its geometrical characteristics are based on the ex-
eter for various fluids[16]. However, the scattering of the
periments of Singh et al.[19] and Choi et al.[20]. A valida-
results confirms the difficulty of describing the evapora-
tion of the model is presented in the last paragraph.
tion accurately. For water, values ranging from 0.01 to 1
have been found. Therefore, an accommodation coefficient
3.1. Standard case equal to 0.1 is chosen for the standard case.
The properties of this “ standard” LHP, having a flat The temperatures at different locations in the LHP are
disk-shaped evaporator, is defined in Table 1. The system plotted in Figure 6 for heat inputs ranging between 10 W
is supposed to be in horizontal orientation. The wick is and 110 W. A temperature difference of several degrees be-
made of nickel and the working fluid is water. The main tween Te and Twe shows the impact of the thermal contact
parameters for the wick are the effective thermal conduc- resistance between the wick and the casing, particularly at
tivity, the accommodation and the contact thermal resis- high heat inputs. Since the mass flow rates in the trans-
tance between the wick and the casing. port lines and the condenser are moderate due to the high
8
 Evaporator design Flat disk-shaped
Wick diameter 40 mm
Wick thickness 3 mm
Wick effective thermal conductivity 5 W · m=1 · K=1
Vapour groove width 1 mm
Wick-wall contact resistance 10=4 K · m2 · W=1
Accommodation coefficient 0.1
Reservoir depth 10 mm
Evaporator diameter 41 mm
Evaporator wall material Stainless steel
Working fluid Water
Fluid charge 7g
Condenser and transport lines i/o diameters 2/2.4 mm
Transport lines length 200 mm
Condenser length 100 mm
Heat sink temperature 22 ◦C
Heat transfer coefficient with the heat sink 2000 W · m=2 · K=1
Ambient temperature 22 ◦C
Heat transfer coefficient with the ambient 5 W · m · K=1
= 2

Table 1: Definition of the standard LHP for the sensitivity analysis

latent heat of vaporization of water, pressure drops in the sink temperature. Above this point, Tc,o starts to increase
loop are not significant. Therefore, the saturation temper- and the shape of all the temperature curves in the LHP is
atures in the reservoir and in the grooves, Tr and Tv , are quasi-linear.
almost equal. The values of Tc,o and Tr,in are also very Figure 7 shows the temperature in the porous structure
similar because the heat transfer between the liquid line in the standard case, for a total heat load Qin of 50 W. For
and the ambient is limited. y = 0 mm, the temperature profile shows a sharp discon-
tinuity at the transition between the groove and the fin.
160 The intense evaporation occurring at x = 0.5 mm causes
Te the convergence of the heatlines toward this point and a
sharp distortion of the lines in its neighbourhood. As y in-
140 Twe
creases, the temperature profile flattens until the temper-
Tv ature is uniform in the x-direction for y > 2 mm. Then,
120 Tr the alternation between the fin and the groove does not
have any influence and the temperature field in the wick
Tc,o in contact with the liquid bulk of the reservoir is uniform.
100
T (°C)

Tr,in As defined in Equation 13, the difference between the

groove temperature and the reservoir temperature is set
80 only by the pressure losses in the transport lines, the
hydrostatic pressure difference and the non-condensable
60 gases partial pressure. Smooth tubes with a relatively
large diameter lead to reduced pressure losses. There-
fore, the LHP operates with a groove temperature and
40 a reservoir temperature almost equal. As a consequence,
the convective cooling due to the liquid flow in the wick is
20 extremely low and heat transfer in the porous structure is
20 40 60 80 100 mainly controlled by heat conduction.
Qin (W) Figure 8 shows the effect of non-condensable gases on
the evaporator temperature. Various quantities of air
Figure 6: Temperature of the LHP in the standard case ranging between 1 µg and 200 µg are simulated in the reser-
voir. For a NCG mass below 10 µg, corresponding to a
A clear transition between the variable conductance partial pressure PNCG equal to about 400 Pa, NCG do not
mode and the fixed conductance mode is observed around have a determinant influence on the evaporator temper-
60 W. For lower heat inputs, the temperature of the liquid ature. However, when the mass of NCG in the LHP is
exiting the condenser is constant and equal to the heat larger, its partial pressure is more important than the sum
9
 Tr rise in the steady-state operating temperature. Besides, it
3 was observed that the performance degrading effect of the
NCG was more pronounced at low heat loads.
79.7 °C 79.7 °C
2.5 3.2. Sensitivity analysis
The LHP operation depends on many distinct param-
2 eters: the geometrical design (size, shape) of the system,
80.2 °C 80.2 °C the thermophysical properties of the working fluid and of
the materials of the loop elements, the heat transfer char-
y (mm)

1.5 acteristics inside and outside the LHP and the pressure
80.7 °C losses in the system. Most of these parameters are easy
80.7 °C to determine if the LHP geometry and the working fluid
1 properties are known. However, some parameters are very
81.2 °C difficult to determine precisely theoretically as it has been
seen in the previous section: the effective thermal con-
0.5 81.7 °C °C
80.7 ductivity of the wick λeff , the contact resistance between
82.2 °C the wick and the evaporator envelope Rc and the accom-
83.2 °C °C
82.7 80.2 °C
0 79.7 °C modation coefficient aev . Furthermore, the heat sink heat
0 Twe 0.5 Tv 1 transfer coefficient hsink and the heat transfer coefficient
x (mm) with the ambient hext are also calculated using correlation
being inherently inaccurate. Thus, in the present section,
a sensitivity analysis is conducted on these parameters to
Figure 7: 2-D temperature field in the wick
see their influence on the model.
The sensitivity coefficient of the function T =
f (xi , xj , xk , ...) in relation to the parameter xi is defined
of the pressure drops in the whole loop. The influence of
as:
NCG becomes then important at low heat input, leading ∂T
to a significant increase of the evaporator temperature. As Si = (72)
∂xi xj ,xk ,...
a consequence, the shape of the characteristic curve of the
LHP is flattened. In order to compare several parameter sensitivities, it is
convenient to define relative sensitivities:
140 ∂T
Si∗ = xi Si = xi (73)
No NCG ∂xi xj ,xk ,...
130 1 µg NCG
10 µg NCG This coefficient enables to quantify the variation ∆T
120 100 µg NCG caused by a relative variation ∆xi /xi of the parameter xi .
The greater the absolute value of the coefficient, the more
200 µg NCG
the function is sensitive to the parameter.
110 400 µg NCG
Unless otherwise mentioned, the sensitivity analysis is
Te (°C)

conducted with the standard LHP defined in Table 1. A

100 slight variation (5 %) of each studied parameter is applied,
the value of the other parameters being constant. This
90 effect on the temperatures at different locations of the LHP
is predicted by the model.
80 Figure 9 presents the relative sensitivity of Te , Tv , Tr,in
and Tc,o to the heat transfer coefficient with the ambient.
70 Sensitivities of Te and Tv are significant for low heat inputs
and decrease with the heat load. Heat transfer with the
60 ambient cools the loop down and leads to a lower opera-
20 40 60 80 100 tional temperature. As expected, hext does not have any
Qin (W) influence on Tr,i and Tc,o in variable conductance mode
because their value is only set by the heat sink. For higher
Figure 8: Influence of the NCGs on the evaporator temperature
heat loads, the increase of the heat losses to the ambi-
ent lead to a lower operational temperature. As a conse-
These results are in accordance with Singh et al.’s exper- quence, a larger part of the condenser is used to conden-
imental study [13]. Their work shows that the net effect of sate the vapour, the subcooling length in the condenser is
the generated NCG in the LHP was to produce an overall smaller and Tc,o increases. The same effect is observed on
10
Tr,in , but partially offset by the heat losses from the liquid 10
line to the ambient.

0
2
-10
0

(dTv/dxi)*xi (K)
-2 -20
(K)

-4 -30
ext
)*h

-6 λeff
ext

-40
(dT/dh

-8 Rc
-50 aev
-10 T
e
hsink
T
v
-12 Tr,in
-60
20 40 60 80 100
-14 Tc,o Qin (W)

20 40 60 80 100 Figure 10: Relative sensitivity of Tv to various parameters

Qin (W)

Figure 9: Relative sensitivity of the LHP temperatures to hext λeff is discussed in the following; it requires a more detailed
analysis on the effect of λeff on the LHP thermal behaviour.
In many experimental configurations, an adequate ther-
mal insulation of the entire LHP enables to reduce con- 20
siderably the heat losses to the ambient. In the following,
this parameter is set to zero in order to better highlight
10
the effect of the other parameters on the LHP.
The relative sensitivity of the vapour temperature to
0
the other parameters is shown in Figure 10. It is clear
that heat transfer inside the evaporator (function of the
(dTe/dxi)*xi (K)

parameters λeff , Rc and aev ) governs the operation in vari- -10

able conductance mode whereas at high heat loads, hsink
becomes the dominant parameter. Indeed, at low heat -20
loads, Tc,o equals Tsink so the heat transfer coefficient with
the heat sink has no influence on the subcooling of the -30
liquid and on the LHP operation in general. In fixed con- λeff
ductance mode, heat transfer in the condenser sets the -40 Rc
LHP operational temperature and the sensitivity of Tv to
hsink is linear. The parameters aev , λeff and Rc have a aev
-50
moderate influence on the vapour temperature and this hsink
influence decreases as the heat load increases. -60
The relative sensitivity of Te to the same parameters is 20 40 60 80 100
shown in Figure 11. The influence of the accommodation Qin (W)
coefficient and of the heat transfer with the heat sink is
almost the same as for Tv . However, the sensitivities to Figure 11: Relative sensitivity of Te to various parameters
Rc and λeff have a different behaviour. While the contact
resistance has a limited effect on Tv , its influence on Te Figure 12 shows the variation of the vapour tempera-
is large and increases with the heat input. Indeed, the ture with the effective thermal conductivity of the wick
difference between Te and Twe is proportional to the heat for several heat loads. An optimal value of λeff is found
transfer rate Qw and to Rc . As the heat input Qin in- between 1 and 2 W · m=1 · K=1 whatever the heat input.
creases, Qw increases almost linearly and leads to a higher This minimum value of Tv is a consequence of the evolution
sensitivity of Te to the contact resistance. of the heat transferred from the casing to the evaporation
The particular shape of the relative sensitivity of Te to zone. As the thermal conductivity of the wick increases,
11
the heat entering the wick becomes larger at the expense 100
of the longitudinal parasitic heat flux. Figure 13 shows the
distribution of heat transfer in the evaporator. The largest
95
part of the heat load Qin enters the wick (Qw ), whereas
the rest is thermally conducted through the evaporator
90
body (Qb ). When the value of λeff is very low, increas-
ing the conductivity enhances both the evaporation (Qev )
85 Q
and the transversal parasitic heat flux (Qw − Qev ). When w

Q/Q (%)
λeff exceeds 1 W · m−1 · K−1 , the increase of the heat flux Qev
≈

in
entering the wick Qw is smaller, leading to an increase of Q -Q
w ev
the transversal parasitic heat flux at the expense of Qev ,
15 Q
that decreases as λeff becomes larger. The maximum evap- b
oration heat flux leads to a minimum vapour temperature
(Figure 12). The same conclusion has been drawn in the 10
numerical study of Siedel et al.[15].
5

150
0
Q=10W 0.5 1 2 3 5 10 20
140 Q=30W λ
-1 -1
(W.m .K )
Q=50W eff

130 Q=70W Figure 13: Distribution of heat transfer in the evaporator (Qin =
Q=90W 50 W)
120 Q=110W

5 W · m=1 · K=1 in the standard case), the sign of the

T (°C)

110 slope of the temperature curve changes as the heat flux

v

increases. As a consequence, the sign of the sensitivity of

100 Te to λeff changes when increasing the heat load, as it is
shown in Figure 11. This is not the case for Tv , since the
90 minimum of the curves Tv (λeff ) does not depend on the
heat flux (Figure 12).
80

80
70
0.5 1 2 3 5 10 20 Q=10W
λ
-1 -1
(W.m .K ) 70 Q=30W
eff Q=50W
Figure 12: Influence of λeff on Tv 60 Q=70W
Q=90W
50 Q=110W
The evolution of the temperature difference Te −Tv with
T -T (K)

λeff for several heat inputs is given in Figure 14. As ex-

pected, this difference decreases when the thermal con-
v

40
ductivity of the wick increases. Indeed, a higher value of
e

λeff leads to lower thermal resistances in the evaporator 30

(Figure 1). Figure 15 presents the evolution of Te with
λeff , which is a consequence of the results obtained in fig- 20
ures 12 and 14. An optimal effective thermal conductivity
exists, for which Te is minimal. However, its value is also 10
dependent on the heat input Qin , contrary to the tem-
perature of the vapour (Figure 12). The optimal value of 0
λeff (obtained for the minimum value of Te ) increases from 0.5 1 2 3 5 10 20
2 W · m=1 · K=1 to 10 W · m=1 · K=1 with the increase of -1 -1
λ (W.m .K )
eff
the heat input from 10 W to 110 W. Therefore, there is
not an optimal value of the thermal conductivity of the Figure 14: Influence of λeff on Te − Tv
wick but a range of optimal values depending on the heat
load. This sensitivity analysis shows that the main parame-
Figure 15 shows that at a given λeff (for example ters of the model are independent. Therefore, their in-
12
 200 between the experimental results and the calculated tem-
Q=10W peratures of the evaporator wall and of the vapour in the
Q=30W grooves. A good agreement is found between the experi-
180 Q=50W mental data and the model. Several unknown parameters
Q=70W were identified by comparison of the predicted values with
Q=90W the experimental data of Singh et al.: the heat transfer co-
160
Q=110W efficient between the condenser wall and the heat sink, the
evaporator wall thickness, the contact resistance and the
Te (°C)

140 accommodation coefficient. As a result, a straight-tube

equivalent condenser is simulated with a heat transfer co-
efficient hsink fixed to 3.2 kW · m=2 · K=1 , considering an
120
outside diameter of 2.4 mm for the tubes. The accommo-
dation coefficient is equal to 0.4. The value of Rc is set to
100 10=5 K · m2 · W=1 . The value of Rc that enables to fit at
best the results is very small, but it has to be noted that Tv
is not really the experimental measurement of the vapour
80 change of temperature, but the temperature of the tube at the exit
slope
of the evaporator. Therefore, these experimental results
0.5 1 2 3 5 10 20 are not sufficient to estimate accurately the parameters.
-1 -1
λ (W.m .K ) Nevertheless, the order of magnitude of the parameters
eff
are close to the one defined for the standard case.
Figure 15: Influence of λeff on Te

100
fluence can be differentiated from each other. Thus, the Te experimental
availability of precise experimental data of several repre- 95
sentative temperatures of the LHP for various heat inputs
Tv experimental
may theoretically lead to a precise determination of these 90 Te model
parameters and provide the model with adequate input
Tv model
parameters. This sensitivity analysis also shows the large 85
influence of the tested parameters on the LHP operation.
Their inaccurate determination can lead to a major error 80
T (°C)

on the LHP operation prediction.

75
3.3. Model validation
The present analytical model is compared to an ex- 70
perimental data set from Singh et al. [19] for a valida-
tion purpose. These authors studied the operational char- 65
acteristics of a flat disk-shaped evaporator LHP, 30 mm
in diameter, using water as working fluid. The 3 mm 60
thick porous wick is made of sintered nickel and its ther-
mal effective conductivity is considered equal to about 55
5.87 W · m=1 · K=1 using Alexander’s formula and based 20 30 40 50 60 70
on Singh et al.’s study[22]: Qin (W)
 −(1−ε)0.59 Figure 16: Comparison between the model and data from Singh et
λl
λeff = λl (74) al. [19]
λwm
where λl and λwm are the thermal conductivities of the
liquid and of the wick material respectively and ε is the
porosity, equal to 75 %. The wick is embedded in a copper 4. Conclusion
evaporator. The vapour and liquid lines, of internal diame-
ter 2 mm, are 150 mm and 290 mm long respectively. A fin- In this paper, a complete analytical model of a LHP has
and-tube condenser, 50 mm long, dissipates heat by forced been developed. Its originality lies in the combination of
convection of air at ambient temperature (i.e. 22 ◦C). energy balance equations for each component of the sys-
No parasitic heat transfer through the evaporator body is tem with analytical solutions for the temperature field in
taken into account, since an O-ring seal prevents heat con- the evaporator. Based on Fourier series expansions, heat
duction to the reservoir. Figure 16 shows the comparison transfer in the wick as well as in the evaporator casing are
13
accurately modelled. This analytical method offers a sim- Rc contact resistance [K.m2 .W −1 ]
ple solution that can be implemented in LHP design anal- S surface area [m2 ]
ysis without the need of large computational resources. Si absolute sensitivity
A sensitivity analysis has been undertaken to assess the Si∗ relative sensitivity
influence of five parameters on the loop operation. This T temperature [K]
analysis enables a better comprehension of the operating T∗ non-dimensional temperature
mechanisms of the LHP as well as a comparative study of V volume [m3 ]
the parameters affecting its temperatures. It appears that X, Y non-dimensional coordinates
these parameters can be experimentally determined using xi,j,k sensitivity parameter
an appropriate test bench with only few temperature mea- x, y axis coordinates [m]
surements.
This model has been validated with a set of experimental Greek Symbols
data from the literature. A good agreement has been met
between the simulation and the experimental results. ∆ difference
The model results show that convection inside the wick ε porosity
does not play a major role and can be neglected. More- λ thermal conductivity [W.m−1 .K −1 ]
over, the temperature field in the wick is almost uniform φ, φ0 heat flux [W.m−2 ]
far from the grooves. ρ density [kg.m−3 ]
In accordance with experimental data from the literature,
the presence of NCGs in the reservoir leads to an increase Subscripts
of the evaporator temperature. This degrading effect is
more pronounced at low heat loads. 2φ two-phase
air air
b evaporator body
Nomenclature
c condenser, convective, contact
A Fourier series coefficient e evaporator
cross-sectional area [m2 ] eff effective
a, a0 , a1 length [m] ev evaporation
aev accommodation coefficient ext external, ambient
B length ratio i inner
Fourier series coefficient in input, inlet
b length [m] l liquid
C Fourier series coefficient NCG non condensable gas
c, c0 , c1 length [m] o outlet, outer
cp specific heat [J.kg −1 .K −1 ] r reservoir
D diameter [m] sen sensible
length ratio sink heat sink
d length [m] sub subcooling
e thickness [m] t total (including convection)
f friction factor v vapour
g gravitational acceleration [m.s−2 ] w wick
H height [m] we wick side of the wick-envelope interface
h heat transfer coefficient [W.m−2 .K −1 ] wm wick material
hlv enthalpy of vaporization [J.kg −1 ]
K coefficient Non Dimensional Numbers
L length [m]
M molar mass [kg.mol−1 ] Bi Biot number
m, n Fourier series increment Pr Prandtl number
mf total fluid charge [kg] Ra Rayleigh number
ṁ mass flow rate [kg.s−1 ] Re Reynolds number
P pressure [P a]
Q heat transfer rate [W ]
q volumetric heat source [W.m−3 ] References

[1] Y. F. Maydanik, Miniature loop heat pipes, in: 13th Interna-

tional Heat Pipe Conference, 2004, pp. 23–35.
[2] Y. Maydanik, Loop heat pipes, Applied Thermal Engineering
25 (2005) 635–657.

14
 [3] T. Kaya, J. Ku, T. T. Hoang, M. K. Cheung, Mathematical
modeling of loop heat pipes, in: 37th AIAA Aerospace Sciences
Meeting and Exhibit, 1999.
[4] A. Delil, V. Baturkin, G. Gorbenko, P. Gakal, V. Ruzaykin,
Modelling of a miniature loop heat pipe with a flat evapora-
tor, in: Proceedings of the 32nd International Conference on
Environmental Systems, 2002.
[5] P. Chuang, An improved steady-state model of loop heat pipes
based on experimental and theoretical analyses, Ph.D. thesis,
Pennsylvania State University (2003).
[6] L. Bai, G. Lin, H. Zhang, D. Wen, Mathematical modeling of
steady-state operation of a loop heat pipe, Applied Thermal
Engineering 29 (2009) 2643–2654.
[7] Y. F. Maydanik, Y. G. Fershtater, N. N. Solodovnik, Loop heat
pipes: Design, investigation, prospects of use in aerospace tech-
nics, Tech. Rep. 941185, SAE International, Warrendale, PA
(Apr. 1994).
[8] S. Launay, V. Sartre, J. Bonjour, Analytical model for char-
acterization of loop heat pipes, Journal of Thermophysics and
Heat Transfer 22 (4) (2008) 623–631.
[9] M. Furukawa, Model-based method of theoretical design anal-
ysis of a loop heat pipe, Journal of Thermophysics and Heat
Transfer 20 (1) (2006) 111–121.
[10] F. Lefèvre, M. Lallemand, Coupled thermal and hydrodynamic
models of flat micro heat pipes for the cooling of multiple elec-
tronic components, International Journal of Heat and Mass
Transfer 49 (7-8) (2006) 1375–1383.
[11] S. Lips, F. Lefèvre, A general analytical model for the design of
conventional heat pipes, to be published, International Journal
of Heat and Mass Transfer.
[12] A. Adoni, A. Ambirajan, V. Jasvanth, D. Kumar, P. Dutta,
Theoretical studies of hard filling in loop heat pipes, Journal of
Thermophysics and Heat Transfer 24 (1) (2010) 173–183.
[13] R. Singh, A. Akbarzadeh, M. Mochizuki, Operational charac-
teristics of the miniature loop heat pipe with non-condensable
gases, International Journal of Heat and Mass Transfer 53 (17)
(2010) 3471–3482.
[14] V. P. Carey, Liquid-vapor phase-change phenomena, Hemi-
sphere, New York, 1992.
[15] B. Siedel, V. Sartre, F. Lefèvre, Numerical investigation of the
thermohydraulic behaviour of a complete loop heat pipe, Ap-
plied Thermal Engineering 61 (2) (2013) 541–553.
[16] I. Eames, N. Marr, H. Sabir, The evaporation coefficient of wa-
ter: a review, International Journal of Heat and Mass Transfer
40 (12) (1997) 2963–2973.
[17] Y. Cao, A. Faghri, Analytical solutions of flow and heat trans-
fer in a porous structure with partial heating and evaporation
on the upper surface, International Journal of Heat and Mass
Transfer 37 (10) (1994) 1525–1533.
[18] F. P. Incropera, D. P. DeWitt, Fundamentals of heat and mass
transfer, John Wiley & Sons, 1996.
[19] R. Singh, A. Akbarzadeh, M. Mochizuki, Operational charac-
teristics of a miniature loop heat pipe with flat evaporator, In-
ternational Journal of Thermal Sciences 47 (2008) 1504–1515.
[20] J. Choi, B. Sung, C. Kim, D.-A. Borca-Tasciuc, Interface engi-
neering to enhance thermal contact conductance of evaporators
in miniature loop heat pipe systems, Applied Thermal Engi-
neering 60 (1–2) (2013) 371–378.
[21] S. Mo, P. Hu, J. Cao, Z. Chen, H. Fan, F. Yu, Effective thermal
conductivity of moist porous sintered nickel material, Interna-
tional Journal of Thermophysics 27 (1) (2006) 304–313.
[22] R. Singh, A. Akbarzadeh, M. Mochizuki, Effect of wick charac-
teristics on the thermal performance of the miniature loop heat
pipe, Journal of Heat Transfer 131 (8) (2009) 082601.
[23] K. Boomsma, D. Poulikakos, On the effective thermal conduc-
tivity of a three-dimensionally structured fluid-saturated metal
foam, International Journal of Heat and Mass Transfer 44 (4)
(2001) 827–836.
[24] C. Li, G. P. Peterson, The effective thermal conductivity of
wire screen, International Journal of Heat and Mass Transfer
49 (21-22) (2006) 4095–4105.

15
Appendix A. Detailed expression of the K-coefficients

 
a0 + a1 λeff 6
K 1 = Sw − + ṁl cp,l 2
(A.1)
2a b πDw
" ∞  
X λeff 1 a0 + a1   mπa 
0
 mπa  
1
−1
−1 
K 2 = Sw 2 2 2
sin mπ cos − cos 1 − e2mπB − 1 − e−2mπB
m=1
a1 − a0 m π 2a a a
   
a0 + a1 λeff a0 + a1 6
+ 1− − ṁl cp,l 2
(A.2)
2a b 2a πDw
" ∞  
X λeff 1 a0 + a1   mπa 
0
 mπa  
1 2mπB −1



−2mπB −1
K3 = −Sw 2 sin mπ cos − cos 1 − e − 1 − e
m=1
a 1 − a 0 m2 π 2 2a a a
 2 #
λeff a0 + a1
− (A.3)
b 2a
  
a0 + a1 λeff 6
K 4 = Sw 1 − − ṁl cp,l 2
(A.4)
2a b πDw
" ∞  
X λeff 1 a0 + a1   mπa 
0
 mπa  
1 2mπB −1



−2mπB −1
K 5 = Sw 2 sin mπ cos − cos 1 − e − 1 − e
m=1
a 1 − a 0 m2 π 2 2a a a
    
a0 + a1 λeff a0 + a1 6
− 1− 1− − ṁl cp,l 2
(A.5)
2a b 2a πDw
" ∞  
X λeff 1 a0 + a1   mπa 
0
 mπa  
1
−1
−1 
K6 = −Sw 2 2 2
sin mπ cos − cos 1 − e2mπB − 1 − e−2mπB
m=1
a1 − a0 m π 2a a a
 
λeff a0 + a1 a0 + a1
+ 1− (A.6)
b 2a 2a
+c1
c0  − c02c +1
  
K7 = hext πc0 c1 −
2 Bib + 1
∞ 
X c 1   mπc 
0
 mπc 
1
− 2λb 2 2
cos − cos
m=1
c1 − c0 m π c c
  mπc   mπc 
1 mπD

πhext c  c   mπc 
0 1
c0 π sin e −1 + −1 + cos + c0 sin
c λb mπ mπ c c
 −1
mπD + Bib
emπD + e−mπD (A.7)
mπD − Bib
  mπc   mπc 
1
πhext c  c   mπc 
0 1
+ −c0 π sin e−mπD − 1 + −1 + cos + c0 sin
c λb mπ mπ c c
 −1 #!
+mπD mπD − Bib
e −mπD
+e
mπD + Bib
 c0 c0 + c1
 1
K8 = hext πc0 c1 −
2 2c Bib + 1
∞ 
X c 1   mπc 
0
 mπc 
1
+ 2λb 2 2
cos − cos
m=1
c 1 − c 0 m π c c
  mπc   mπc 
1
πhext c  c   mπc 
0 1
c0 π sin emπD − 1 + −1 + cos + c0 sin
c λb mπ mπ c c
 −1 !
mπD −mπD mπD + Bib
e +e (A.8)
mπD − Bib
  mπc   mπc 
1
πhext c  c   mπc 
0 1
+ c0 π sin e −mπD
−1 + −1 + cos + c0 sin
c λb mπ mπ c c
16
  −1 #!
+mπD mπD − Bib
e −mπD
+e
mπD + Bib
 
 c0  φ0 d c0 φin 1
K9 = hext πc0 c1 − − Text
2 λb cφ0 Bib + 1
∞ 
X 1  mπc 
0
+ 2φin c 2 2 sin
m=1
m π c
  mπc   mπc 
1
πhext c  c   mπc 
0 1
−c0 π sin emπD − 1 − −1 + cos + c0 sin
c λb mπ mπ c c

−1 
1 + e2mπD (A.9)
  mπc  
1
πhext c  c   mπc 
0
 mπc 
1
+ −c0 π sin e−mπD − 1 + −1 + cos + c0 sin
c λb mπ mπ c c

−1 i
1 + e−2mπD
+c1
c0  − c02c +1
  
K10 = hext πc0 c1 −
2 Bib + 1
∞ 
X c 1   mπc 
0
 mπc  πh
1 ext c
 c   mπc 
0
 mπc 
1
− 2λb 2 π2
cos − cos −1 + cos + c 0 sin
m=1
c 1 − c 0 m c c λ b mπ mπ c c
" −1  −1 #!
mπD −mπD mπD + Bib +mπD mπD − Bib
× e +e + e −mπD
+e (A.10)
mπD − Bib mπD + Bib
 c0  c0 + c1 1
K11 = hext πc0 c1 −
2 2c Bib + 1
∞ 
X c 1   mπc 
0
 mπc  πh
1 ext c
 c   mπc 
0
 mπc 
1
+ 2λb 2 π2
cos − cos −1 + cos + c 0 sin
m=1
c 1 − c 0 m c c λ b mπ mπ c c
" −1  −1 #!
mπD −mπD mπD + Bib +mπD mπD − Bib
× e +e + e −mπD
+e (A.11)
mπD − Bib mπD + Bib
 
 c0  φ0 d c0 φin 1
K12 = hext πc0 c1 − − Text
2 λb cφ0 Bib + 1
∞ 
X 1  mπc0 πhext c  c 
  mπc 
0
 mπc 
1
− 2φin c 2 2 sin −1 + cos + c0 sin (A.12)
m=1
m π c λb mπ mπ c c
h
−1
−1 i
× 1 + e2mπD − 1 + e−2mπD
   
π (Lc − L2φ ) πLl
K13 = exp − exp − (A.13)
ṁl cp,l (1/hl Dc,i + 1/hsink Dc,o ) ṁl cp,l (1/hl Dl,i + 1/hextDl,o )
   
πLl π (Lc − L2φ )
K14 = exp − 1 − exp − Tsink
ṁl cp,l (1/hl Dl,i + 1/hextDl,o ) ṁl cp,l (1/hlDc,i + 1/hsinkDc,o )
  
πLl
+ 1 − exp − Text (A.14)
ṁl cp,l (1/hl Dl,i + 1/hextDl,o )

Appendix B. Second order equation for the determination of Tv

0 =Tv2 [cp,l (1 − K13 ) (K6 (K1 + K2 + K7 ) − (K4 + K5 ) (K3 + K8 ) + K8 Rc /Sc (K6 (K1 + K2 ) − K3 (K4 + K5 )))]
+ Tv [cp,l (1 − K13 ) ((K3 + K8 ) K4 ∆T + (K3 K4 − K1 K6 ) K8 ∆T Rc /Sc + K6 (−K1 ∆T − K7 ∆T + K9 − Qin ))
+ (hlv − cp,l K14 ) (K6 (K1 + K2 + K7 + K8 (K1 + K2 )) − (K4 + K5 ) (K3 + K8 + K3 K8 Rc /Sc )) (B.1)
+hlv (K11 (K1 + K2 + K7 + K3 K7 Rc /Sc ) − (K3 + K8 + K3 K8 Rc /Sc ) (K10 + hext Sr ))]
+ cp,l K14 (K6 (Qin − K9 ) + ∆T (K6 (K1 + K7 + K1 K8 Rc /Sc ) − K4 (K3 + K8 + K3 K8 Rc /Sc )))
17
+ hlv ∆T ((K3 + K8 ) (K4 + K10 ) − (K6 + K11 ) (K1 + K7 ) + K3 K8 (K4 + K10 ) Rc /Sc − (K3 K7 K11 + K1 K6 K8 ) Rc /Sc )
+ hlv ((K3 + K8 + K3 K8 Rc /Sc ) (Qin − K12 + hext Sr (∆T + Text )) − (K6 + K11 + K3 K11 Rc /Sc ) (Qin − K9 ))

18