1. INTRODUCTION
Heat pipes have remarkable heat transfer capacity, high thermal conductivity, and are easy to maintain.
These devices can extract heat and transport it to a condensation zone by changing the fluid phase. They are also highly efficient and their design is simple [1]. Heat pipes are composed of three zones: evaporator, adiabatic zone, and condenser.
Such pipes have a wick through which a working fluid flows due to capillary pressure. A wickless heat pipe is called a thermosyphon, in which the working fluid flows due to gravity. Therefore, thermosyphons must operate vertically or at low inclinations.
Heat transfer in a thermosyphon is caused by the evaporation-condensation of the working fluid. The device has liquid in the evaporation zone; the heat input in this zone evaporates the liquid, and the steam goes up to the condensation zone through the adiabatic section. In the condenser, heat is deflected to the condensation fluid that remains at low temperatures. Then, the liquid film from the condensation returns to the evaporator through the pipe walls due to gravity [1].
Due to their high efficiency and simple operation, heat pipes are used in many engineering applications [2], such as heating, air conditioning [3], electronic devices [4], humidity control [5], geothermal applications, and others.
Therefore, it is essential to analyze the behavior of thermosyphons and their efficiency under different operating conditions.
It is difficult to experimentally describe the changing phase in a thermosyphon. For that reason, it is recommendable to study the behavior of thermosyphons using computational fluid dynamics (CFD) software. Numerical simulations achieve a huge accuracy and can reduce the time and cost of a study [6].
Some researchers have studied the effect of heat input, fill ratio, inclination angle, working fluid, and how these parameters affect the performance of two phase closed thermosyphons (TPCTs).
Regarding the fill ratio, Alizadehdakhel et al. [1] observed that a thermosyphon’s efficiency increases at fill ratios from 0.3 to 0.5. At fill ratios above 0.5, the vapor and liquid film flows increase, but the device thermal resistance increases, too. They obtained a maximum efficiency in the TPCT with a heat input of 500 W and a fill ratio of 0.5. They observed that, with heat inputs higher than 500 W, the liquid film thickness rises, increasing the thermal resistance.
Pulsed heat input has been studied by Kafeel et al. [7], who examined the effect of increasing the heat input between 10 % and 20 %. They observed that the system stabilizes 200 s after every increase. With low increases of around 10 %, they found that the operating temperature rises around 14 K; in turn, with 20 % increases, they found rises of 11 K.
Fadhl et al. analyzed the temperature distribution in a TPCT [2]. They simulated a TPCT with water as working fluid and studied the temperature in different device zones. They found that higher heat inputs in the evaporation section produce an increase in the temperature distribution and the convection coefficient in the evaporation zone. Also, they found a ratio between the effective thermal conductivity and the efficiency of the device, and that the thermal resistance increases when the heat input is low.
Another important phenomenon in TPCTs is the geyser effect, which occurs when the liquid fill ratio is high and the heat input is low. Jouhara et al. [8] observed that when the geyser effect occurs the temperature has a cyclic behavior due to big steam bubbles that drag fluid above them to the condensation zone. Those steam bubbles were observed in their CFD simulation.
Regarding the working fluid, other authors have analyzed distilled or deionized water, nanofluids, and refrigerants. Fadhl et al. [9] studied the effects of using two refrigerants: R134a and R404a.
They found that these fluids present the same behavior in a CFD simulation, and the vapor bubbles are smaller than those formed with distilled water as working fluid. Therefore, the thermal resistance of the device decreases. Ong et al. [10] experimentally studied the behavior of the R410a refrigerant as working fluid. The experimental procedure was carried out with low heat inputs (between 20 W and 100 W), increasing 20 W every 30 minutes. They used water was used as condensate fluid with a constant mass flow of 0.25 kg/s. The inclination angles varied between 30 º and 90º; and the fill ratio, between 0.5 and 1.
The effects of inclination and the fill ratio were considerable at heat inputs greater than 60 W. The highest heat transfer coefficient in condensation zone was 1200 W/m2K with a fill ratio and heat input of 0.75 and 100 W, respectively. The lowest thermal resistance was 0.17 K/W, produced with a heat input of 100 W and a fill ratio of 1.
TPCTs can operate with low angle inclinations. Zhanget et al. [6] studied the effects of angle inclination and working fluid wettability on thermosyphons. When the contact angle is low, the temperature in the wall pipe is low too, and the bubbles move away from de wall quickly, thus increasing the heat transfer. As a result, when the angle contact is larger than 90 °, the steam bubbles adhere to the wall longer, thus increasing the thermal resistance. They found that increasing the inclination angle reduces the temperature measurement fluctuation.
Additionally, some geometric configurations of TPCTs can be studied using CFD simulations. Fertahi et al. [11] conducted a numerical simulation of a TPCT using a smooth condenser and a condenser with fins. In their simulation, they observed that, using the condenser with fins, the heat transfer coefficient was higher than with its smooth counterpart because the vapor was accumulated in the fins and remained longer, thus increasing the heat transfer to the condenser fluid.
Wang et al. used CFD simulations to analyze the effects of the volumetric flow rate of the condensate fluid and the heat input on a radial wickless heat pipe [12].
Through CFD simulations by Wang et al. [12]. In their simulation, the liquid film from the phase change was observed at around 20 s; and the higher the volumetric flow rate of the condensed fluid, the higher the heat transfer of the device. Also, their CFD simulation showed an average thermal resistance between 0.024 K/W and 0.033 K/W under stable operating conditions. Furthermore, nucleate boiling and condensation film were the dominant transfer mechanism in the radial heat pipes they simulated. Wang et al. [13] studied the effect of ammonia as working fluid in a Loop heat pipe (LHP) made of carbon steel. Their analysis used both CFD simulation in ANSYS Fluent 16.0 and experiments. In the simulation, the boundary condition in the evaporator was a constant temperature, and the operating conditions were studied at 13 °C, 18 °C, and 22 °C in the evaporation zone, and between 2 °C and 10 °C using water as condensing fluid. The simulations showed that the heat transfer is the highest with a 200 ml/min volumetric flow rate of condensing fluid and the lowest temperature of the condensing fluid.
The effect of the condensing fluid on the temperature is significant at volumetric flow rates up to 100 ml/min.
Controllable loop heat pipes (CLHPs) are thermosyphons with control valves in the vapor and liquid lines. These devices are used in refrigerators to extract heat from their compartments. The effect of non- condensable gases on CLHPs has been studied by Cao et al. [14] using CFD simulations. They used R134a refrigerant as working fluid, air as non-condensable gas, and a propylene glycol solution as condensing fluid. They conducted their analysis to understand the operating behavior of CLHPs with non-condensable gases. Their results showed that the operating conditions with 0.47 % air as non-condensable gas are acceptable. Their CFD simulation showed that the heat transfer in the device decreases from 380.6 W to 254.5 W when the proportion of air rises from 0 to 0.62 %, and the heat transfer dropped by 72 % with 0.62 % non-condensable gases. The effect of non-condensable gases was negligible up to 0.21 % of the mass of air in the device.
Azzolin et al. analyzed a TCTP with integrated water storage through numerical simulations [15]. Their mathematical model was analyzed with MATLAB Simulink, and the results were validated with CFD simulations and experimentally. Two different diameters were analyzed, 8 mm and 10 mm; the inclination angle was between 15º and 90º; and the input heat flux was constant, 800 W/m2. The simulation showed that the temperature at the top of the tank was the highest with an 8-mm diameter pipe at all the inclination angles, and the temperature in the middle was the lowest with an inclination of 15º. At the bottom of the tank, the temperature was the highest with the 10-mm diameter pipe and inclination angles between 30º and 90º. The highest efficiency was 67.2 %, produced with the 10 mm diameter pipe and an inclination angle of 45º.
Zhang et al. studied an air conditioning system for data centers using simulations [16] that implemented a distributed-parameter method and several operating conditions: constant temperature in the evaporation zone, 27 ºC; outdoor temperature, between 12 ºC and 22 ºC; and air volumetric flow rate, between 0.35 m3/s and 1.25 m3/s. They found that the heat transfer from the device to the air is proportional to the air volumetric flow rate the temperature difference between the evaporation and condensation zones. Also, their simulations showed that the effects of diameter and length of the evaporator are negligible for the heat transfer if the temperature difference remains constant.
Wang et al. implemented the Lee model to simulate heat and mass transfer in CFD software [17]. They also modified such model to simulate the geyser effect in a TPCT. Those modifications were implemented in ANSYS Fluent 16.0 through a User-Defined-Function (UDF).
The implementation of the UDF in the model aimed to properly simulate the nucleation phenomenon according to the overheating temperatures with the modification of the relaxation factor in the equations. In the numerical simulation with the original Lee model, the liquid pool temperature in the evaporation zone remained close to the saturation temperature, while in the numerical simulation with the modified model the liquid pool and the steam bubbles presented overheating temperatures. The simulation with the modified model represented more accurately the bubble formation observed experimentally.
Despite the number of numerical simulations performed to represent the operation of a TPCT, the literature lacks studies about the most suitable models and parameters to reproduce the different heat and mass transfers occurring in this type of devices because these phenomena are hard to observe experimentally. Therefore, this study aims to deepen our understanding of some models available in a commercial CFD code (ANSYS Fluent), which could be applied to the numerical simulation of a TPCT in order to obtain an accurate representation of the physical phenomena governing the operation of such device.
Particularly, the models and parameters of turbulence, density, phase change, and phase interfaces were modified and the results were analyzed to determine the feasibility of the application of those models and parameters to the simulation of a TPCT.
2. METHODOLOGY
2.1 Theoretical model
ANSYS Fluent software was used for the numerical simulation of a TPCT. The Volume of Fluid (VOF) model was used to simulate the multi-phase flow. All the phases present in the simulation were considered immiscible [18]. The VOF model is based on the fact that the sum of the fraction volumes of all the phases present in a computational cell equals one [1] (1):
Where ∝𝑙 and ∝𝑣 are the liquid fraction and vapor volume, respectively. If the cell is full of liquid, then ∝𝑙= 1; but, if it’s full of steam, ∝𝑣= 1.
The density in a computational cell is defined as (2):
Where 𝜌𝑣 and 𝜌𝑙 correspond to the vapor and liquid density, respectively. The continuity equation for VOF model is (3), (4):
Where 𝑢⃗ represents the average velocity, and 𝑚𝑣 and 𝑚𝑙 are the mas flows due to evaporation-condensation respectively. The momentum equation including the volume fraction terms is set as (5):
Where I represent the unit tensor, 𝑝 denotes the pressure and 𝑔 the gravity 𝜇 denotes the viscosity, this value is in function of volume fraction in a computational cell as (6):
FCSF determine the volumetric surface tension This term was proposed by Brackbill [19] and it has the following expression (7):
Where 𝜎𝑣𝑙 is the surface tension coefficient (8).
The energy equation is:
Where 𝑇 is the mixture temperature. The thermal conductivity in each cell is calculated through the volume fraction (9):
Where 𝑘𝑙 and 𝑘𝑣 are the energy transfer rates of liquid and steam, respectively. The specific energy of the system is (10):
The terms 𝑒𝑙 and 𝑒𝑣 are calculated in function of the specific heat (11) and (12):
Where 𝑐𝑝,𝑙 and 𝑐𝑝,𝑣 are the specific heats for liquid and vapor, respectively.
Fluent´s mass and energy transfer model is the Lee model [20]. The model has equations for calculate the mass and energy source terms from each phase as follows (13), (14), (15), (16):
The energy transfer between phases is (17) and (18):
Where ∆𝐻 is the vaporization enthalpy.
2.2 Development of the simulation
The geometry of the TPCT, the mesh and the simulation was performed in ANSYS 18.1. The model geometry was drawn in DesingModeler. A 2D model was proposed and the specifications are shown in Table 1.
The mesh was designed with Meshing in ANSYS Workbench. A cut cell method was used to structure the mesh. The cells have an initial spacing of 0.1 mm and a growth factor of 1.2 near the inner wall to simulate the condensate bubbles.
The resultant mesh is shown in Fig. 1.
A total of 58,800 grids were used for the computational calculation. The wall thickness was created as a boundary condition in ANSYS Fluent, as shown in Fig. 1.
In Fluent, a transient solution scheme was selected to simulate the operation of the TPCT. The pressure-based solver was selected for this simulation. The VOF model was chosen for the numerical simulation with two Eulerian phases.
A temperature-dependent density correlation was proposed according to data in liquid and steam tables based on values obtained from [21], (19):
The vapor density was calculated using the ideal gas law. The surface tension coefficient was defined through (20) based on data obtained from [22], (20):
2.2.1 Boundary conditions
A constant heat flux of 6973.925 W/m2 was imposed in the evaporation zone. In the top and the bottom caps, the imposed heat flux was 0, i.e., they were considered adiabatic surfaces. The condenser zone wall The condenser zone wall had a constant output heat flux of 9176.22 W/m2.
A non-slip boundary was defined on the inner wall of all the zones in the TPCT.
2.2.2. Solution methods
The SIMPLE algorithm was chosen for the Pressure-Velocity coupling, and the schemes for the discretization were first-order upwind for momentum and energy equations, Geo- Reconstruct for volume of fraction equation, and PRESTO for pressure. The time step size was 5x10-4 s, and the convergence criteria was based on residuals around 10-6 or lower for energy and 10-3 or lower for the other equations.
2.2.3 Initial parameters
Some initial parameters and models were ideal gas density, a laminar viscosity model, a mass and energy transfer coefficient of 0.1, and sharp interface modelling. However, changes were made to these parameters and models during the simulation to achieve convergence and fitting of the experimental data.
For example, viscosity was shifted to a urbulent k-Ω model, sharp interface modelling was shifted to disperse interface, and the mass and energy transfer coefficient was varied between 0.1 and 100. Other parameters and models remained unchanged.
2.3 Experimental procedure
To validate the CFD simulation, an experimental evaluation was conducted.
The experimental set up is shown in Fig. 2. type J. Five thermocouples were located in the evaporation zone, three thermocouples in the condensation zone, and two thermocouples in the inlet and outlet of the condensate fluid.
A pressure transducer was used in the device to monitor the pressure, and a submersible pump condenser increased, too. The effect of this difference on the thermal performance of the TPCP was considered using the thermal resistance, which was calculated employing the average temperatures in the device’s zone at each supplied power through (21):
Was used to ensure a constant volumetric flow rate of cooling fluid in the condensation zone. Also, thermostatic bath was used to ensure a constant temperature of 25 °C for the cooling fluid.
The power input in the evaporator was simulated through a direct current power supply and an electrical resistance.
The power input was calculated using Watt’s law, the power input was calculated using Watt’s law and the measurements of current and voltage.
The output power in the condensation zone was calculated using Newton’s heat transfer law.
3 RESULTS
3.1 Experimental evaluation
Four tests were conducted with the same conditions: an initial vacuum pressure of 71.11 kPa, an input power of 10W initially, and a 10 W increase every 35 minutes up to 70 W.
The results in Fig. 3 show that the temperature rises as the input power increases. After each power increase, the temperature rose up to a steady state after 19 minutes. Between 19 and 35 minutes, the changes in average temperature in the device’s zones were not significant, between 19 and 35 minutes, the changes in average temperature in the device’s zones were not significant, and it is considered the saturation temperature because an increase in the pressure inside the device due to evaporation caused an increase in the saturation temperature. In Fig. 3, the temperature curve on the bottom shows the average temperatures in the condensation zone, while the upper curve represents the average temperatures in the evaporation zone.
The experimental results in Fig. 4 show that the thermal resistance of the TPCT decreased as the power input increased. Furthermore, the thermal resistance shows an asymptotic behavior close to 0.5 (°C/W). In Fig. 3, the temperature difference between the evaporation and condensation zones is greater at high input powers. Equation (19) shows a proportionality with such difference, while input power is inversely proportional to thermal resistance. Then, as can be seen in Fig. 4, the higher input power was the most influential factor in the increase in temperature differences in the device.
3.2 CFD model validation
For the validation, the thermosyphon in the CFD simulation had the geometry shown in Table 1. 50 W were imposed as the boundary condition in the evaporation and condensation zones. The time step was 5x10-4 s.
Along the length of the thermosyphon, there was a good agreement in temperature values between the experiments and the CFD simulations.
The results are shown in Fig. 5. The deviation of 25 °C in the central zone is believed to be due to the fact that the numerical model is not able to predict the geyser phenomenon that occurs in the real thermosyphon, which helps to stabilize the wall temperature. Therefore, future work should aim to overcome this limitation of numerical simulations.
3.3 Numerical Simulation
A numerical simulation with phase change is a fluid flow problem of special attention since when the immiscible phases are dependent each other, it is easy to obtain divergence. That is why it is necessary to properly define the physical phenomenon there. A series of numerical simulations were performed to determine the behavior of the fluid flow, the average temperatures and the heat transfer in the device. The process consisted in the variation of some parameters mentioned below, while the other parameters remained constant. The parameter varied where:
3.3.1 Turbulence model
ANSYS Fluent has different models such as k-Ɛ model and k-ɷ model for simulation of turbulence. Two regimes were analyzed through the laminar and kɷ models, for laminar and turbulent regime, respectively. The results obtained for both models were the similar. Besides, Fluent allows to observe the kinetic energy dissipation. For the turbulent model implemented, the kinetic energy dissipation has a low value as can be seen in the contour in the evaporation zone in Fig. 6. The highest value 8.82x10-3 m2 /s2 was obtained in the liquid pool, while, in other zones in the device, the turbulent kinetic energy is negligible. In Fig. 7, it can be seen the velocity contours for the laminar and turbulence viscous model, where it is possible to observe that the velocity profile with turbulence viscous model corresponds to low velocities between 0 and 10-1 meters per second, which confirms the velocity scale of the laminar viscosity model.
3.3.2 Density model
Density plays an important role when the buoyancy is present inside the physic phenomenon. In a TPCT, the saturated vapor goes up in the center of the pipe while the liquid from condensation returns to the bottom of the device on the tube walls. Due to high temperatures, it is possible that the vapor reaches high velocities, causing an increase of Courant´s number. This number is the ratio between the computational cell length and the flow velocity. When this number is greater than 0.25, the simulation diverges. By the other hand, high speeds produce, in the same way, high pressures, causing divergence in equations or poor convergence.
In the present study, two simulations were performed with two density models for vapor: constant density and ideal gas law density. With constant density, the vapor velocity rising is low, but the high temperatures predominate in evaporation section where there is no liquid. Also, with constant density model for vapor, a maximum speed of 0.7 meters per second is observed in the evaporation zone, while, in the condensation zone, the vapor reaches low speeds, which prevents the temperature homogenization in the condensation zone and increases the temperature differences between the device´s zones, which increases the thermal resistance.
Using the ideal gas law for the gas phase, the temperature is evenly distributed along the pipe, and the velocity of vapor is higher than with constant steam density. Fig. 8a and Fig. 8b, shows the temperature contours for both cases.
3.3.3 Energy and mass transfer coefficient calculated using the phase change
The 0.1 energy and mass transfer coefficient used in the Lee model [13]-[18] controls the mass flow in the evaporation-condensation phenomenon and the energy source in each phase. With high coefficients for mass flow and energy source in the condensation zone, the vapor temperatures in zones adjacent to the wall pipe and the center of the pipe may remain near the saturation temperature.
Imposing a constant temperature as boundary condition in condensation zone, it is possible to obtain a constant output heat rate, and the cool temperature remains constant only in the cell near the wall pipe.
With coefficients between 0.1 and 1, the simulation remains in a steady state, but it is not possible to see a condensate film.
Around the first seconds, the output heat rate is high due to temperature differences.
Then, the vapor reaches thermal equilibrium with the wall, causing a low output heat rate. Therefore, the coefficient of mass flow and energy sources for evaporation must always be low because, otherwise, it causes divergence.
Fig. 9 shows the temperature profile produced with different mass and energy transfer coefficients in phase change.
3.3.4 Interface Modelling
Two options for interface modelling were evaluated: sharp and disperse interface. Choosing the right type of interface is important because this parameter controls the mixing of phases and the movement of each phase relative to the other. In a disperse interface, two or more phases are interpenetrating, while the sharp interface is appropriate when a distinct interface is present between the phases. Fig. 10 and Fig. 11 were obtained using the disperse and sharp interfaces, respectively. In the case of the disperse interface, vapor was mixed with the liquid phase in all the zones of the device when the phase change occurred, while the sharp interface restricted the condensing fluid to the walls of the device. Therefore, in order to reproduce the phase change phenomenon in a TPCT, a sharp interface is highly recommended.
4. CONCLUSIONS
A numerical simulation of the operation of a TPCP was performed to visualize the effect of choosing some parameters offered in ANSYS Fluent.
According to the results, using a turbulence viscous model is not necessary because the velocity profiles were the same with laminar and turbulent models.
Additionally, using a turbulence model requires a lot of computation time.
Therefore, the physical phenomenon inside the TPCT can be better represented with a laminar viscous model. In terms of the density model, a variable density model improved the temperature distribution inside the device because it allowed for vapor buoyancy.
The phase change during condensation should be an isothermal process for vapor that remains at the center of the pipe.
With high mass and energy transfer coefficients during condensation, the vapor remains close the saturation temperature. Finally, the disperse interphase model did not show good results because the phases were mixed when the phase changed, while the sharp interphase model allowed the phases to be immiscible with a different interface between both phases, which is closer to the physical phenomenon that takes place inside TPCTs.