View metadata, citation and similar papers at core.ac.

uk brought to you by CORE

provided by Brunel University Research Archive

Numerical Modelling of the Temperature Distribution in a Two-

Phase Closed Thermosyphon

Bandar Fadhl 1,2, Luiz C. Wrobel 1, Hussam Jouhara 1

1
School of Engineering and Design, Brunel University, Uxbridge, Middlesex UB8 3PH, UK
2
On leave from Department of Mechanical Engineering, Umm Al-Qura University, Makkah, KSA

Keywords: Thermosyphon, Computational fluid dynamics (CFD), Multiphase flow, Phase change
material, Evaporation, Condensation

ABSTRACT
Interest in the use of heat pipe technology for heat recovery and energy saving in a vast
range of engineering applications has been on the rise in recent years. Heat pipes are playing a
more important role in many industrial applications, particularly in improving the thermal
performance of heat exchangers and increasing energy savings in applications with commercial
use. In this paper, a comprehensive CFD modelling was built to simulate the details of the two-
phase flow and heat transfer phenomena during the operation of a wickless heat pipe or
thermosyphon, that otherwise could not be visualised by empirical or experimental work.
Water was used as the working fluid. The volume of the fluid (VOF) model in ANSYS
FLUENT was used for the simulation. The evaporation, condensation and phase change
processes in a thermosyphon were dealt with by adding a user-defined function (UDF) to the
FLUENT code. The simulation results were compared with experimental measurements at the
same condition. The simulation was successful in reproducing the heat and mass transfer
processes in a thermosyphon. Good agreement was observed between CFD predicted
temperature profiles and experimental temperature data.

1 INTRODUCTION

A heat pipe is a two-phase heat transfer device with a highly effective heat transfer rate
through evaporating and condensing a fluid that is circulating in a sealed container. A wickless
heat pipe, or a two-phase closed thermosyphon, relies on gravitational forces to return the
working fluid to the evaporator. This is different from a wicked heat pipe, where the working
fluid is returned from the condenser by capillary forces [1-4]. Heat pipes have been
successfully used for waste heat energy recovery in a vast range of engineering applications,
such as heating, ventilation, and air conditioning (HVAC) systems [5], ground source heat
pumps [6], water heating systems [7] and electronics thermal management [8]. This is mainly
because of their simple structure, special flexibility, high efficiency, good compactness, and
excellent reversibility [9-12]. Thermosyphons have three sections, which are the evaporator at
the bottom end, where heat is added and the liquid is vaporised, the condenser at the top end,

1
where heat is released and the vapour is condensed, and an adiabatic section in the middle
between the evaporator and condenser [13].
In a thermosyphon, heat is added to the evaporator where a liquid pool exists, changing the
liquid into vapour. The high temperature and pressure cause the vapour to flow and pass
through the adiabatic section toward the condenser. The vapour adjacent to the condenser’s
wall gives up its latent heat that is absorbed in the evaporator section. The condensed liquid is
then transported back to the evaporator due to gravity [14].
Two-phase closed thermosyphons have been extensively used in many applications [15].
However, only a limited number of CFD numerical simulation studies on two-phase closed
thermosyphons have been published. Alizadehdakhel et al. [1] provided a two-dimensional
model and experimental studies in which they investigated the effect of input heat flow and fill
ratio of the working fluid on the performance of a two-phase closed thermosyphon. They
validated their study using experimental results. Three input heat flow rates of 700, 500, and
350W and three fill ratios of 0.3, 0.5, and 0.8 were considered. Under these operating
conditions, they found the performance of the thermosyphon improved when the input heat
flow was increased from 350 to 500W. Further, they discovered the best performance was at a
fill ratio of 0.5. The authors reported a term called “heat performance”, which they calculated
by using the following equation for different fill ratios:

Qout
η= × 100
Qin

However, this term is not usual in heat pipe publications. In general, the thermal performance
term used to characterize thermodynamics at different heat throughputs is the total thermal
resistance.
Legierski et al. [14] provided CFD modelling and experimental measurements of heat
and mass transfer in a horizontal wicked heat pipe. They investigated the effectiveness of the
heat pipe thermal conductivity in a transient state during start-up of the pipe operation and
during temperature increases. The authors used a heat pipe that was 200 mm long with 4 mm
diameter and 25 mm length for the evaporator and condenser. They also used two containers,
one for hot water (90 oC) at the evaporator section and one for cold water (ambient
temperature) at the condenser section. They developed a three-dimensional CFD model to
simulate the internal vapour flow. They found that the effective thermal conductivity of the
wicked heat pipe depended on the time in the range between 15 x 103 and 30 x 103 W/m K, and
achieved its steady-state value after approximately 20 to 30s. However, the authors did not
consider in the CFD modelling the phase change material from liquid phase to vapour phase, as
well as condensation in the condenser section and pool boiling in the evaporator section.
Zhang et al. [16] developed a two-dimensional heat and mass transfer model for a disk-
shaped flat two-phase thermosyphon used in electronics cooling. The authors simulated the
vapour flow inside the flat two-phase thermosyphon as a single-phase flow. They compared
their predicted model with experimental results to determine the factors that affected the axial

2
thermal resistance of a thermosyphon. This model was limited as it considered the flow inside
the flat thermosyphon as a pure vapour phase only.
Joudi and Al-Tabbakh [17] numerically studied a two-phase thermosyphon solar
domestic hot water system, by using a computer simulation. They used R-11 as a working fluid
in the thermosyphon. Firstly, the authors validated the computer program and calculation
procedure by comparing the results with those obtained with single-phase systems. They then
performed calculations for the two-phase thermosyphon system. In their calculations, they
evaluated mass flow rate, saturation pressure, and temperature in the collector and condenser,
together with tank temperature and collector and condenser thermal efficiencies. The results of
the study showed that the collector efficiency of the two-phase system was approximately 20%
greater than in a single-phase system. Further, the response of the two-phase system in
reaching maximum tank temperature and efficiency was faster than in a single phase system.
This study was only a mathematical model and did not include any flow visualisation.
Annamalai and Ramalingam [18] carried out an experimental investigation and CFD
analysis of a wicked heat pipe using ANSYS CFX. The authors considered the region inside
the heat pipe as a single phase of vapour and a wick region as the liquid phase. They compared
the predicted surface temperature along the evaporator and condenser walls and the vapour
temperature inside the heat pipe with the experimental data. This model treated the flow inside
the heat pipe as a single-phase and did not include the evaporation, condensation and phase
change processes.
De Schepper et al. [19] developed a model to simulate the evaporation process of a
hydrocarbon feedstock in a heat exchanger. They used the VOF and UDF techniques to
simulate flow boiling including the phase change process. They proposed correlations to
calculate the mass and heat transfer between the phases that were able to simulate the
evaporation and boiling phenomena inside the convection section of a steam cracker. This
model was for the convection section in a steam-cracking furnace; however, it did not include
the heat pipe system.
Lin et al. [20] built a CFD model to predict the heat transfer capability of miniature
oscillating heat pipes. The effects of different heat transfer lengths and inner diameters at
different heat inputs were used to analyse the heat transfer capability of MOHPs. They
compared the predicted model with experimental results. This model did not visualise the
internal phenomena of evaporation, condensation and phase change inside the MOHPs.
Heat pipe technology is currently still under development. However, there are limited
studies on the validation of predictions for modelling closed two-phase thermosyphons or
wickless heat pipes. Further, a CFD simulation of a wickless heat pipe that considers all the
details of heat transfer phenomena inside the heat pipe has not yet been reported. Hence, a gap
still exists for further CFD work to model a wickless heat pipe. Additionally, CFD models can
reduce the amount of experimental work. Therefore, in this paper, a comprehensive CFD
modelling has been employed to cover all details of two-phase flow and heat transfer
phenomena during the operation of a straight wickless heat pipe. Moreover, a user-defined
function (UDF) has been used to complete the FLUENT code in order to simulate the phase
change material.

3
2 EXPERIMENTAL APPARATUS

In order to validate the CFD findings, an experimental apparatus was built to carry out a
thermal performance investigation on a typical wickless heat pipe.
The experimental apparatus used in the current investigation is shown in Figure 1. The
apparatus consists of a two-phase closed thermosyphon (TPCT), a rope heater, the cooling
water circuit, and instrumentation. The apparatus was fixed on a framework to insure vertical
orientation under all test conditions.
The TPCT was manufactured from a 22mm outer diameter, 0.5m-long smooth copper
tube with a wall thickness of 0.9mm. It consists of a 0.2m-long evaporator section, a 0.1m-long
adiabatic section and a 0.2m-long condenser section (see Figure 1).
The evaporator section was heated by a rope heater with a maximum power output of
500W at 220V, which was evenly wrapped and not directly positioned above any of the
thermocouples that are used to measure the surface temperature of this section. The energy
output of the heater was controlled by a variac. The evaporator section was wrapped in a layer
of fire-proof insulation before it was wrapped with suitable thermal insulation layers to
minimise any heat losses to the ambient. The condenser section was cooled using a double pipe
concentric heat exchanger with an insulated outer surface. The 10-cm long adiabatic section
was also well insulated to ensure no heat energy interactions take place with the ambient. The
insulated adiabatic section wall temperature was used as an indicator of the TPCT working
temperature. The TPCT was charged with triple-distilled and degassed water at a 50%
evaporator filling ratio.
The cooling water circuit provided the heat exchanger of the condenser section with
water at predefined conditions for the inlet temperature and flow velocity. The water was
supplied through a constant-head water circuit to ensure constant mass flow rate through the
condenser’s heat exchanger. The heat from the water leaving the heat exchanger was removed
in a secondary reservoir, using a chilled water coil, before it was returned to the main tank
using a dedicated pump. A proportional-integral-derivative control system (PID) was used to
control the chiller system to ensure water supply at the heat exchanger inlet was constant
throughout the testing programme. A flow meter and a valve arrangement were used to control
and measure the inlet volume flow rate into the shell of the heat exchanger.
The experimental apparatus was equipped with calibrated instrumentations to measure the
power throughputs, temperatures and flow rate data.
The temperature distribution along the TPCT was monitored using eight thermocouples.
As shown in Figure 1, two thermocouples, labelled as Te1 and Te2 were used to monitor the
evaporator section and were placed 40mm and 160mm from the bottom. Another
thermocouple, labelled as Ta, was positioned at the centre of the adiabatic section, while the
condenser section was monitored using five, evenly spaced, thermocouples, labelled as Tc1 to
Tc5. These five thermocouples were also used to confirm the non-existence of non-condensable

4
gases (NCGs) within the heat pipes throughout the conducted tests. Two additional
thermocouples were also used to monitor the input and the output water temperatures from the
heat exchanger. These two thermocouples were stainless steel shielded and were positioned at
the centre of the flow using two compression fittings.
All of the thermocouples were K-type (NiCr/NiAl) and were read and monitored using a
32-channel DataScan system, which was connected to a dedicated PC for an online data
recording at 1Hz scanning frequency. The cooling water flow rate was measured using an
inline flow meter.

3 EXPERIMENTAL PROCEDURE

At the start of each experiment, and after fixing the water flow rate to the desired value,
the electrical heat input was set at its minimum level of 50W. The equipment was then allowed
to stabilise for 30 minutes prior to any readings being taken. Temperature readings, from all the
rig thermocouples, were then monitored for two hours, using the datalogger. The two-hour
monitoring period was designed to ensure that no degrading of performance was taking place
during operation. The procedure was repeated for various electrical input powers between 50
and 500W with an increment of 50W. Each test was repeated three times for each power
setting to confirm repeatability. The repeated tests were done after the whole power range was
covered to confirm stable thermal characteristics of the TPCT after a prolonged operational
time.

3.1 Data reduction

The effective overall thermal resistance of the TPCT was calculated by applying the
electrical analogue in the form:

∆Te −c
REXP = (1)
Q

T e1 +T e 2
∑ Tcj
j = 1→5
∆Te −c = Te _ av − Tc _ av Te _ av = Tc _ av =
here ; 2 and 5 are the average internal wall
temperatures in the evaporator and condenser sections, respectively, and Q is the power
throughput. The internal wall temperatures of the TPCT were measured after considering the
thermal conduction across the shell wall.

5
3.2 Uncertainty analysis

The main source of uncertainty for the calculated R came from the temperature readings,
which were measured using K-type thermocouples with a measurement uncertainty of
±(0.05% rdg + 0.3°C) and 0.5% rdg for the power readings.
According to Taylor [21], the propagation of uncertainties associated with the calculated R
( S R ) can be estimated from the equation:

2 2
 S ∆T   SQ 
=S R R  e −c  +   (2)
 ∆Te −c   Q 

where:

=
S ∆Te−c STe _ av 2 + STc _ av 2 : The uncertainty associated with ∆Te −c

SQ in : The uncertainty associated with the reading of the energy throughput Q

By calculating SR, for the entire experimental range, the maximum uncertainty associated
with the resulting R values was found to be around 5.5%, which is an acceptable value in
engineering applications.

4 MODEL DESCRIPTION

In this model, the commercial code ANSYS FLUENT 13.0 and the Volume of Fluid
(VOF) method have been applied for the modelling of a closed two-phase thermosyphon.
There are two main approaches to the numerical calculations of multiphase flows, which are
the Euler-Lagrange approach and the Euler-Euler approach. The Euler-Lagrange approach
treats the fluid phase as a continuum and a dispersed phase such as bubbles or droplets as a
second phase, in which the dispersed-phase volume fractions are not to exceed 10%. As the
current application considers the volume fraction of the second phase exceeds 10%, the Euler-
Euler approach has been adopted as it uses the idea of phasic volume fraction in which the
volume of a phase cannot be occupied by the other phases. These volume fractions are assumed
to be continuous functions of space and time [22].

4.1 Volume of fluid (VOF) model

Numerical solutions based on the finite volume method are more difficult for multiphase
flows than for a single-phase flow. The reasons for this difficulty are that the interfaces
between the phases are not stationary and physical properties such as density and viscosity

6
change at the interfaces between the different phases, which requires an intensive
computational effort. The volume of fluid (VOF) technique, therefore, has been used to solve
these problems by determining the motion of all phases and defining the motion of the
interfaces indirectly from this result [19,23,24].
The VOF technique can be applied to model two immiscible fluids with a clearly defined
interface between the phases, and is used for surface-tracking applied to a fixed mesh. In the
VOF model, one set of Navier-Stokes equations are solved through the computational domain
and used to track the motion of the different phases by defining the volume fraction of each
phase [22]. The VOF model relies on the fact that each cell in the domain is occupied by one
phase or a combination of the two phases. In other words, if αL is a volume fraction of liquid
and αV is a volume fraction of vapour, the following three conditions are possible:

• αL = 1 : The cell is fully occupied by liquid

• αL = 0 : The cell is fully occupied by vapour
• 0 < αL < 1 : The cell is at the interface between the liquid and vapour phases

When the third condition occurs, the volume fractions of all phases sum to unity [22].

4.2 Navier-Stokes equations for VOF model

The governing equations of mass continuity, momentum and energy are used to describe
the motion of the working fluid in a thermosyphon. This will be explained in the next section.

4.2.1 Continuity equation for VOF model (Volume fraction equation)

By applying the physical principle of conservation of mass to the fluid, the continuity
equation has the following form:

 ∂ρ
∇. (ρu ) = − (3)
∂t
where ρ is the density, u is the velocity and t is the time.
Solution of the above equation for the volume fraction of one of the phases is used to
track the interface between the phases. Thus, the continuity equation of the VOF model for the
secondary phase (L) can be expressed as:

 ∂
∇. (α L ρ L u ) = − (α L ρ L ) + S m (4)
∂t

where Sm is the mass source term used to calculate the mass transfer during evaporation and
condensation.
The continuity equation shown above can be called the volume fraction equation and
this equation will not be solved for the primary phase; the primary-phase volume fraction is
computed based on the following constraint:

7
 n

∑α
L =1
L =1 (5)

When the cell is not fully occupied by the primary phase (V) or the secondary phase
(L), a mixture of the phases L and V exist. Thus, the density of the mixture is given as the
volume-fraction-averaged density and takes the following form:

ρ = α L ρ L + (1 − α L )ρV (6)

4.2.2 Momentum equation for VOF model

The forces acting in the fluid were considered to be gravitational, pressure, friction and
surface tension. In order to consider the effect of surface tension along the interface between
the two phases, the continuum surface force (CSF) model proposed by Brackbill et al. [25] has
been added to the momentum equation

α L ρ L CV ∇α V + α V ρV C L ∇α L
FCSF = 2σ LV (7)
ρ L + ρV

where σ LV is the surface tension coefficient and C is the surface curvature.

By taking into account the above forces, the momentum equation for the VOF model
takes the following form:

∂ 
∂t
( )
(ρu ) + ∇. (ρu u ) = ρg − ∇p + ∇. µ ∇u + ∇u T − 2 µ∇. uI  + FCSF (8)
 3 

where g is the acceleration of gravity, p is the pressure, and I is the unit tensor.
The momentum equation depends on the volume fraction of all phases through the
physical properties of density and viscosity. Thus, the dynamic viscosity µ is given by

µ = α L µ L + (1 − α L )µV (9)

A single momentum equation is solved throughout the computational domain, and the
calculated velocity is shared among the phases.

4.2.3 Energy equation for VOF model

The energy equation for the VOF model has the following form:

8
 ∂
(ρe ) + ∇. (ρeu ) = ∇. (k . ∇T ) + ∇. ( pu ) + S E (10)
∂t

where SE is the energy source term used to calculate the heat transfer during evaporation and
condensation.
The VOF model treats the temperature T as a mass-averaged variable and the thermal
conductivity is calculated as:

k = α L k L + (1 − α L ) kV (11)

The VOF model also treats the internal energy e as a mass-averaged variable in the following
form:

α L ρ L e L + α V ρV eV
e= (12)
α L ρ L + α V ρV

where eL and eV are based on the specific heat cp of the phase and the shared temperature, given
by the caloric equation of state:

e L = c p , L (T − Tsat ) (13)

eV = c p ,V (T − Tsat ) (14)

A single energy equation is also solved throughout the domain for both phases, and the
calculated temperature is shared among the phases.

4.3 Mass and heat transfer during the evaporation and condensation processes

FLUENT does not have the ability to simulate the phase change material during the
evaporation and condensation processes. In order to circumvent this problem, a user-defined
function (UDF) has been used to complete the existing FLUENT code. This UDF is essentially
required to calculate the mass and heat transfer between the liquid and vapour phases during
the evaporation and condensation processes, determined by the source terms in the governing
equations, particularly the continuity and energy equations. Source terms proposed by De
Schepper et al. [19] have been used to calculate the mass and energy transfer. Mass sources, SM
in the volume fraction equation and energy sources, SE in the energy equation used in the
present work can be found in Table 1, where Tmix and Tsat are the mixture temperature and
saturation temperature, respectively, and LH stands for latent heat.

9
 Mass and energy sources in Table 1 have been implemented in the UDF and linked to the
governing equations in FLUENT. The volume fraction for each phase in the cell has been
defined by the VOF model. Therefore, the evaporation process required two mass sources for
the calculation of the mass transfer, Eq. (15) describing the amount of mass taken from the
liquid phase and Eq. (16) describing the amount of mass added to the vapour phase. The same
procedure takes place for the condensation process, Eq. (17) and Eq. (18) describing the
amount of mass transfer from vapour to liquid phase.
For heat transfer, a single source term for both phases is required in the evaporation or
condensation. Calculation of heat transfer has been determined by multiplying the mass source
with the latent heat for evaporation or condensation, Eq. (19) and Eq. (20), respectively.
Furthermore, it can be seen in Table 1 that the temperature is introduced as a mixture
temperature rather than liquid or vapour temperature. The reason as mentioned before is that
the VOF model associates some variables such as temperature and velocity with the mixture
phase, not with a specific phase.

5 SIMULATION MODEL

A two-dimensional model was developed to simulate the two-phase flow and heat
transfer phenomena in a thermosyphon. A total length of 500 mm of a closed copper tube, as
can be seen in Figure 1, is used as the thermosyphon geometry, with 22 and 20.2 mm for the
outer and inner diameters, respectively. According to the experimental condition, the
thermosyphon was divided into three sections represented by the evaporator and condenser
sections, with an adiabatic section between them. Both evaporator and condenser have 200 mm
length, while the adiabatic section has 100 mm length.
The temperature distribution along the outer wall of the thermosyphon was monitored
using eight different positions, which are the thermocouple positions as shown in Figure 2.
According to the experimental setup, Te1 and Te2 were used to record the average temperature
of the evaporator section, while Tc1 to Tc5 were used to record the average temperature of the
condenser section. Ta was used to record the average temperature of the adiabatic section.
Different mesh sizes were used to test grid independence. The average temperature of the
evaporator (Tevap av) and condenser (Tcond av) sections for different mesh sizes were monitored
and are shown in Table 2. For the heating power of 172.87 W, it was found that almost the
same temperature differences between the evaporator and condenser sections were obtained for
different mesh sizes. As a result, the mesh size of 69,092 Quad, Map cells was selected for the
simulation analysis. Near the left and right walls, five layers of cells are used in order to
capture the thin liquid film that develops near the wall, as shown in Figure 3. One cell layer has
been used for the upper and bottom walls, as no heat conduction is considered through these
walls.

10
6 BOUNDARY CONDITIONS

A non-slip boundary condition is imposed at the inner walls of the thermosyphon. In

order to simulate the heating and evaporation, a constant heat flux is defined at the wall
boundaries of the evaporator section, depending on the power input. A zero heat flux is defined
as boundary condition on the adiabatic section, assuming this section is insulated. The
condenser section is cooled as a result of heat released when vapour condenses. It is assumed
that the condenser is cooled by water, according to the experimental apparatus. Thus, a
convection heat transfer coefficient is defined as boundary condition on the condenser’s wall.
The corresponding heat transfer coefficients have been calculated using the formula:

Qc
hc =
2πrLc (Tc ,av − T∞ )
(21)

where hc is the condenser heat transfer coefficient, Qc is the rate of heat transfer from the
condenser, Lc is the condenser height, Tc,av is the condenser average temperature and T∞ is the
average temperature of the condenser cooling water. The values of T∞, Qc and Tc,av in the
above equation and in Table 3 come from the experiments. The condenser heat transfer
coefficients are determined based on the experimental data, and are shown in Table 3.
The effect of surface tension along the interface between the two phases is included by
using the following equation, driven from the steam table:

σ LV = 0.09805856 − 1.845 × 10 −5 T − 2.3 × 10 −7 T 2 (22)

where T is the shared temperature.

The model considered water as the working fluid with a fill ratio of 0.5 (the ratio of
initial liquid volume per total volume of the evaporator) and the following equation driven
from the steam table is used for their density:

ρ L = 859.0083 + 1.252209T − 0.0026429T 2 (23)

Figure 4 illustrates the boundary conditions implemented to the computational model.

7 SOLUTION STRATEGY AND CONVERGENCE CRITERION

A transient simulation with a time step of 0.0005 s is carried out to model the dynamic
behaviour of the two-phase flow. The time step has been selected based on the Courant
number, which is the ratio of the time step to the time a fluid takes to move across a cell. For
VOF models, the maximum Courant number allowed near the interface is 250 [26]. For a time

11
step of 0.0005, the Courant number is less than 3. The simulation reaches a steady state after
around 60 s.
A combination of the SIMPLE algorithm for pressure-velocity coupling and a first-order
upwind scheme for the determination of momentum and energy is included in the model. Geo-
Reconstruct and PRESTO discretization for the volume fraction and pressure interpolation
scheme, respectively, are also performed in the simulation. In the present studies, the numerical
computation is considered to have converged when the scaled residual of the mass and velocity
components is less than 10-4.
Water vapour is defined as the primary (vapour) phase and water liquid is defined as the
secondary (liquid) phase. For the calculation of the mass and heat transfer during the
evaporation and condensation processes, a temperature of 373K is used as the boiling
temperature and the latent heat in the UDF code is 2455 kJ/kg. When the simulation is started,
the liquid pool in the evaporator is heated first. Once the saturation temperature (373K) is
reached, evaporation starts and phase change occurs. The saturated vapour is then transported
upward to the condenser, where it condenses along the colds walls forming a thin liquid film.

8 SIMULATION RESULTS

The simulation results of the boiling and condensation processes in the thermosyphon
reached a quasi-steady state after around 60 s. The temperature distribution along the outer
wall of the thermosyphon for heating powers of 172.87 W and 376.14 W are shown in Tables 4
and 5, respectively. Eight different positions have been used to monitor the average
temperature for the evaporator, adiabatic and condenser sections. Table 6 shows the surface
average temperatures in the evaporator (Teav), adiabatic (Ta) and condenser (Tcav) sections, in
addition to the thermal resistance of the system and the relative error (RE) between CFD
simulation and experimental results (EXP). The simulation results of the VOF model showed
the same trend as the experimental data. The average relative error of evaporator, adiabatic and
condenser average temperatures are 7.9%, 9.9% and 1.9%, respectively.
Figure 5 shows the experimental and simulation results of the outer surface temperature
distribution along the thermosyphon for different heat inputs. The distance between 0 and 200
mm indicates the evaporator section. The distance between 200 and 300 mm indicates the
adiabatic section, while the distance between 300 and 500 mm indicates the condenser section.
The predicted CFD evaporator average temperature has deviated from the experimental results
due to the consideration of a continuous heat power input along the length of the evaporator
section where, in the experiment, a wire heater is evenly wrapped around the evaporator
section to ensure it was not directly above a thermocouple. As shown in Figure 5, the
condenser section temperature shows better agreement with the experimental results. As a
result of no heat loss in the adiabatic section, the temperature is raised in the surface of this
section due to the axial conduction heat transfer.

8.1 Performance of a closed two-phase thermosyphon

The performance of a thermosyphon can be characterised by the overall thermal

resistance. The overall rate of heat transfer to the system 𝑄𝑄̇ is proportional to the effective

12
temperature difference between the heat source to the evaporator and the heat sink from the
condenser, and inversely proportional to the equivalent thermal resistance to heat transfer
between the two regions. The overall rate of heat transfer can be defined as:

∆T
Q = (25)
R

Hence, the effective overall thermal resistance of a thermosyphon R is calculated using the
following equation:

Te av CFD − Tc av CFD
RCFD = (26)
Qin

where Teav CFD and Tcav CFD are the average temperature in the evaporator and condenser,
respectively, and Qin is the heating power input. Different heating power inputs have been used
to investigate the performance of a closed two-phase thermosyphon. Figure 6 illustrates that
the predicted thermal resistance is in good agreement with the experimental data as the thermal
resistance of the thermosyphon decreases with increasing heating power load. For heating
powers above 170 W, the thermal resistance stays relatively independent of the heating power
input rate. For lower heating inputs, the thermal resistance tends to increase. In summary, the
CFD simulation results are able to show the variation trends of the thermal performance of the
thermosyphon for different heat throughputs.

9 FLOW VISUALISATION OF CFD SIMULATION

Figures 7 and 8 show the volume fraction contours of pool boiling in the evaporator and
condensed liquid film in the condenser, respectively, for a heating power of 172.87 W. A red
colour illustrates the presence of only vapour (vapour volume fraction = 1), while a blue colour
stands for the presence of only liquid (vapour volume fraction = 0). In Figure 8, focus is made
on the condensed liquid film in the lower region of the condenser. At the beginning of the
process, the liquid pool that initially filled half of the evaporator was heated by imposing a
constant heat input. At positions where the liquid reached the boiling temperature, the liquid
starts to evaporate and phase change occurs as shown in Figure 7. This continuous evaporation
of liquid results in a decrease in the liquid volume fraction and an increase of the vapour
volume fraction. At those positions where the liquid evaporates, bubbles are formed and
transported toward the top region of the liquid pool.
Following the above process, saturated vapour is transported upward to the condenser.
As the vapour reaches the condenser’s wall, where a convection heat transfer coefficient
boundary condition is defined, as shown in Table 3, the vapour condenses along the cold walls
forming filmwise condensation as shown in Figure 8. This liquid will then fall back to the
evaporator section and recharge the liquid pool.
As shown in Figure 9, the temperature contours during the thermosyphon operation
have been recorded. The heating power is 172.87 W. At the beginning, the temperature in the
evaporator section increased due to the heating power input, as shown in Figure 9 (0.5 s and 1

13
s). When the temperature of the evaporator section reached the boiling temperature, the phase
change begins. The region of high temperature in the evaporator section expands due the
vapour moving upward, as shown in Figure 9 (1.5 s to 3 s). As the heating power in the
evaporator section continues, the vapour flows across the adiabatic section to the condenser
section, as shown in Figure 9 (4 s and 5 s). Then, a high temperature region appears in the
condenser section due to the vapour reaching this section, as shown in Figure 9 (10 s). The
high temperature of the condenser section starts to decrease, corresponding to vapour
condensing to liquid and, with the help of gravity, the condensed liquid falls back to the
evaporator section. The above cycle describes the process of heat transfer during the operation
of the thermosyphon. After that, the temperature distribution inside the thermosyphon becomes
uniform as shown in Figure 9 (30 s to 60 s).

10 CONCLUSIONS

The main objective of this work is the development of a CFD model that allows to perform
simulations of the evaporation and condensation phenomena in a thermosyphon. The
simulation of these processses is one of the steps required to model the complete system in
order to consider the phase change material by implementing the appropriate source terms in
the flow governing equations. These source terms, determining the mass and heat transfer
between the liquid and vapour phases, have been linked to the main hydrodynamic equations of
FLUENT.
The CFD results of this work show that FLUENT with the VOF method can successfully
model the complex phenomena inside the thermosyphon. From the flow visualisation, it is
found that the CFD simulation was able to reproduce the operation of the thermosyphon,
including the pool boiling in the evaporator section and the condensed liquid film in the
condenser section.
The average surface temperature along the thermosyphon has been compared with the
experimental results at the same condition, showing that the predicted results agreed with the
experimental results quite well. The thermal performance of the thermosyphon has also been
characterised at different heat throughputs by the effective overall thermal resistance, and it is
found that increasing the heating power inputs above 172 W has improved the thermal
performance of thermosyphon.

Acknowledgment
The first author deeply appreciates the financial support provided by the Saudi Cultural
Bureau in London, the Ministry of Higher Education and the Mechanical Engineering
Department, Umm Al-Qura University.

14
References
[1] Alizadehdakhel A, Rahimi M, Alsairafi AA. CFD modeling of flow and heat transfer in a
thermosyphon. International Communications in Heat and Mass Transfer 2010; 37: 312-318.

[2] ESDU. Heat pipes - general information on their use, operation and design. ESDU Manual 80013
1980.

[3] Faghri A. Heat Pipe Science and Technology, Taylor & Francis: Washington, D.C., 1995.

[4] Dunn P, Reay DA. Heat Pipes, Pergamon Press: New York, 1994.

[5] Kerrigan K, Jouhara H, O'Donnell GE, Robinson AJ. Heat pipe-based radiator for low grade
geothermal energy conversion in domestic space heating. Simulation Modelling Practice and Theory
2011; 19: 1154-1163.

[6] Jouhara H, Meskimmon R. Experimental investigation of wraparound loop heat pipe heat exchanger
used in energy efficient air handling units. Energy 2010; 35: 4592-4599.

[7] Mathioulakis E, Belessiotis V. A new heat-pipe type solar domestic hot water system. Solar Energy
2002; 72: 13-20.

[8] Weng Y-, Cho H-, Chang C-, Chen S-. Heat pipe with PCM for electronic cooling. Applied Energy
2011; 88: 1825-1833.

[9] Jouhara H. Economic assessment of the benefits of wraparound heat pipes in ventilation processes
for hot and humid climates. International Journal of Low-Carbon Technologies 2009; 4: 52-60.

[10] Parand R, Rashidian B, Ataei A, Shakiby K. Modeling the transient response of the thermosyphon
heat pipes. Journal of Applied Sciences 2009; 9: 1531-1537.

[11] Ochsner K. Carbon dioxide heat pipe in conjunction with a ground source heat pump (GSHP).
Applied Thermal Engineering 2008; 28: 2077-2082.

[12] Du J, Bansal P, Huang B. Simulation model of a greenhouse with a heat-pipe heating system.
Applied Energy 2012; 93: 268-276.

[13] Cengel YA. Heat Transfer : A Practical Approach. McGraw-Hill: Boston; Toronto, 2003.

[14] Legierski J, Wiecek B, de Mey G. Measurements and simulations of transient characteristics of

heat pipes. Microelectronics and Reliability 2006; 46: 109-115.

[15] Jiao B, Qiu LM, Zhang XB, Zhang Y. Investigation on the effect of filling ratio on the steady-state
heat transfer performance of a vertical two-phase closed thermosyphon. Applied Thermal Engineering
2008; 28: 1417-1426.

[16] Zhang M, Liu Z, Ma G, Cheng S. Numerical simulation and experimental verification of a flat two-
phase thermosyphon. Energy Conversion and Management 2009; 50: 1095-1100.

15
[17] Joudi KA, Al-Tabbakh AA. Computer simulation of a two phase thermosyphon solar domestic hot
water heating system. Energy Conversion and Management 1999; 40: 775-793.

[18] Annamalai AS, Ramalingam V. Experimental investigation and computational fluid dynamics
analysis of a air cooled condenser heat pipe. Thermal Science 2011; 15: 759-772.

[19] De Schepper SCK, Heynderickx GJ, Marin GB. Modeling the evaporation of a hydrocarbon
feedstock in the convection section of a steam cracker. Computers & Chemical Engineering 2009; 33:
122-132.

[20] Lin Z, Wang S, Shirakashi R, Winston Zhang L. Simulation of a miniature oscillating heat pipe in
bottom heating mode using CFD with unsteady modeling. International Journal of Heat and Mass
Transfer 2013; 57: 642-656.

[21] Taylor JR. An Introduction to Error Analysis: The Study of Uncertainties in Physical
Measurements. University Science Books: Sausalito, 1997.

[22] ANSYS FLUENT Theory Guide (Release 13.0). Multiphase Flows. ANSYS, Inc., November
2010, (chapter 17), pp. 455-568.

[23] Versteeg HK, Malalasekera W. An Introduction to Computational Fluid Dynamics; The Finite
Volume Method. Prentice-Hall: Harlow, Second ed., 2007.

[24] Anderson JD. Computational Fluid Dynamics The Basics with Applications. McGraw-Hill: New
York, 1995.

[25] Brackbill JU. A continuum method for modeling surface tension. Journal of Computational
Physics 1992; 100: 335-354.

[26] ANSYS FLUENT User Guide (Release 13.0). Modelling Multiphase Flows. ANSYS, Inc.,
November 2010, (chapter 26), pp. 1143-1144.

16
 Tables

Table 1: Construction of mass and energy sources [19]

Thermal Phase Change Temp.
Phase Source Term
Energy process Condition
Tmix − Tsat
Liquid S M = −0.1ρ La L (15)
Tsat
Evaporation Tmix > Tsat
T − Tsat
Vapour S M = 0.1ρ La L mix (16)
Mass Tsat
Transfer T − Tmix
Liquid S M = 0.1ρV aV sat (17)
Tsat
Condensation Tmix < Tsat
T − Tmix
Vapour S M = −0.1ρV aV sat (18)
Tsat
Tmix − Tsat
Evaporation Tmix > Tsat S E = −0.1ρ La L LH (19)
Heat Tsat
Transfer Tsat − Tmix
Condensation Tmix < Tsat S E = 0.1ρV aV LH (20)
Tsat

Table 2: Grid-independence results

Mesh size (cells) 19,603 69,092 87,800
Tevap av (K) 378.71 378.37 378.19
Tcond av (K) 326.54 326.96 327.79
R CFD(K/W) 0.3017 0.2974 0.2915

Table 3: Condenser heat transfer coefficients for different heat inputs

Condenser
Evaporator cooling Condenser
water
Qin T∞ Qc Tc av hc
W K W K W/m2.K
100.41 298.9 95.1 312.41 509.3
172.87 301.45 162.6 318.07 707.6
225.25 302.95 192.2 320.55 790.1
275.60 305.2 236.6 325.95 824.9
299.52 306.3 254.8 323.91 1046.6
376.14 309.4 336.6 330.33 1163.5

17
Table 4: Comparison between experimental data and CFD simulation for heat input of 172.87 W
TEXP TCFD RE TEXP av TCFD av RE av
Section Position
K K % K K %
Te1 345.75 378.33 9.42
Evaporator 341.6 378.37 10.78
Te2 337.45 378.40 12.14
Adiabatic Ta 327.45 362.41 10.68 10.68
Tc1 320.55 329.54 2.80
Tc2 318.85 326.54 2.41
Condenser Tc3 317.95 325.95 2.52 318.07 326.96 2.80
Tc4 317.05 325.64 2.71
Tc5 315.95 327.13 3.54

Table 5: Comparison between experimental data and CFD simulation for heat input of 376.14 W
TEXP TCFD RE TEXP av TCFD av RE av
Section Position
K K % K K %
Te1 376.75 385.14 2.23
Evaporator 370.2 385.05 4.01
Te2 363.65 384.97 5.86
Adiabatic Ta 342.75 370.11 7.98 7.98
Tc1 328.95 327.12 0.56
Tc2 325.55 323.66 0.58
Condenser Tc3 332.45 323.15 2.80 330.33 323.96 1.92
Tc4 331.35 322.70 2.61
Tc5 333.35 323.17 3.05

Table 6: Comparison between experimental data and CFD simulation for different heat inputs
Thermal
Source Evaporator Adiabatic Condenser
Resistance
Teav Teav Ta Ta Tcav Tcav
Qin RE RE RE REXP RCFD
EXP CFD EXP CFD EXP CFD
W K K % K K % K K % K/W K/W
100.41 343 376.18 9.67 321.25 363.25 13.07 312.412 328.35 5.10 0.3046 0.4763
172.87 341.6 378.37 10.76 327.45 362.41 10.68 318.07 326.96 2.80 0.1361 0.2974
225.25 348.1 379.92 9.14 331.05 364.94 10.24 320.55 323.47 0.91 0.1223 0.2506
275.6 356.1 381.6 7.16 335.55 365.62 8.96 325.95 327.36 0.43 0.1094 0.1967
299.52 358.75 382.41 6.60 336.25 365.46 8.69 323.91 324.81 0.28 0.1163 0.1923
376.14 370.2 385.06 4.01 342.75 370.11 7.98 330.33 323.96 1.93 0.1060 0.1624
Average relative error % 7.89 9.94 1.91

18
 Figures

Figure 1: The experimental apparatus

19
 Tc5

Tc4

Tc3

Tc2

Tc1

Ta

Te2

Te1

Figure 2: Geometry and

dimensions

20
y

x
19,603 69,092 87,800

Figure 3: A section of the computational mesh

21
 y

x
Figure 4: Boundary conditions of the thermosyphon

22
 400
Qin = 100.41 W 400
380 Qin = 172.87 W Qin = 100.41
380
Temperature (K)
Temperature (K)

Qin = 225.25 W W
360 Qin = 275.6 W 360 Qin = 172.87
Qin = 299.52 W
340 W
340 Qin = 376.14 W
320
320
300
300 0 100 200 300 400 500 600
0 100 200 300 400 500 600 Length (mm)
Length (mm) (b) Simulation
(a) Experiment

Figure 5: Temperature comparison between experiments and simulations for different heat inputs

23
 0.6
CFD
Thermal resistance (K/W)

0.5 EXP
0.4

0.3

0.2

0.1

0.0
100.41 172.87 225.25 275.6 299.52 376.14
Heating power (W)

Figure 6: Relationship between thermal resistance and heating power

24
 Evaporator
y

x
Time(s)
Figure 7: Contours of volume fraction of pool boiling in the
evaporator section at different times

25
Lower region of condenser
Adiabatic

Time(s)
Figure 8: Contours of volume fraction of condensed liquid film in the lower region of
condenser at different times

26 y

x
Temp. (K)

x
Time(s)

Figure 9: Temperature contours at different times

27