An experiment and three-dimensional numerical simulation of pulsating heat pipes

International Journal of Heat and Mass Transfer 150 (2020) 119317

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/hmt

An experiment and three-dimensional numerical simulation of pulsating heat pipes Duy-Tan Vo, Hyoung-Tak Kim, Junghyuk Ko, Kwang-Hyun Bang∗ Department of Mechanical Engineering, Korea Maritime and Ocean University, 727 Taejongro, Yeongdogu, Busan 49112, Republic of Korea

a r t i c l e

i n f o

Article history: Received 18 August 2019 Revised 11 December 2019 Accepted 4 January 2020

Keywords: Pulsating heat pipe Two-phase flow CFD Fluent

a b s t r a c t An experimental and analytical study of pulsating heat pipe (PHP) has been carried out in order to develop a more reliable numerical simulation model. The test PHP was made of transparent Pyrex tubes with the inner diameter of 1.85mm and a total of sixteen tubes formed eight turns of PHP. Both heating and cooling were provided by water jackets in which inlet and outlet temperatures were precisely measured for estimating heat transfer rate. The working fluid was R123. Visualization using a high-speed camera showed various flow patterns as well as the fluid motions. The wall temperature measured at various locations revealed its relation to the fluid motion and the direction of circulation. The circulation motion was dominant in most tests. The heat transfer rate was measured and the difference for filling ratio of 50 and 60% was little. A three-dimensional computational fluid dynamics modeling has been developed for pulsating heat pipe using ANSYS Fluent. The VOF model with variable density and vapor pressure relation successfully simulated the circulating motions of PHP. The predicted wall temperatures showed the same indication of mode of flow motion and the flow direction as observed in the experiment. The predicted heat transfer rates agreed well with the experimental data within 5%. The success of the simulation of the experiment implies that the realizable k−ε turbulence model is appropriate to PHP simulation and the use of variable density for liquid and vapor and the vapor pressure equation is crucial. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Heat pipe is, in general, a device through which heat transfers better with the help of latent heat of vaporization in such that the working fluid charged inside the device evaporates in the heating section and the vapor moves to the cooling section and it condenses. In conventional heat pipes, the condensed liquid returns to the heating section with the help of capillary force enhanced with so-called wicks fabricated on the inside surface of heat pipe. Therefore, the wick design governs the rate of returning condensate flow and thus this limits the heat transfer performance of the conventional heat pipes. Another type of heat pipes is an oscillating heat pipe, or sometimes called a Pulsating Heat Pipe (PHP). It has no wick, but plain inner surface. The pulsating heat pipe is an efficient heat transport device that employs evaporation of charged working fluid at the heat source and condensation at the heat sink of the heat pipe. The pressure difference between the evaporation and condensation parts acts to move the working fluid inside the heat pipe in an os∗

Corresponding author. E-mail address: [email protected] (K.-H. Bang).

https://doi.org/10.1016/j.ijheatmasstransfer.2020.119317 0017-9310/© 2020 Elsevier Ltd. All rights reserved.

cillatory pattern. A construction of multiple turns out of a long, capillary-size tube can amplify this pressure difference. While the heat transport performance of the conventional heat pipe is limited mainly by the capillary force in narrow flow area through wicks, in the PHP, there is always a driving force of pressure difference through relatively large flow area of the whole tube cross section by different saturation pressures corresponding to heating and cooling temperature as long as the heating part and the cooling part are kept at different temperature. Despite the advantages of simple, higher thermal performance and economical construction of a pulsating heat pipe, a practical application to cooling of electronic packages has been scarce. This is mainly because the pulsating action may not sustain in some poor designs or in abnormal conditions. The pulsating action can slow down or seize sometimes. With the rapid development of semiconductor technology, the operating power of micro-electro processors such as in smartphones and mobile devices has continued to increase higher and higher. This trend has brought a high demand of an efficient heat dissipation technology that must prevent hot spots. The pulsating heat pipe with high performance is a promising heat transfer device for such a demand.

2

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

Nomenclatures cp k m˙ p q T v

specific heat, J/kg K thermal conductivity, W/mK mass flow rate, kg/s pressure, Pa heat flux, W/m2 temparature, K velocity, m/s

Greek symbols α void fraction μ viscosity, Ns/m2 ρ density, kg/m3 Subscript in l out rad sat v

inlet liquid outlet thermal radiation saturation vapor

In the past several decades, the experimental as well as theoretical studies on the pulsating heat pipe have been abundant since Akachi [1] patented it in 1990. The physical parameters that influence the operation and the performance of the pulsating heat pipes are known to be tube diameter, thermal conditions (i.e., temperature range), type of working fluid, charging ratio, number of channels, and the arrangement (or inclination against gravity direction). Clear understanding and quantitative evaluation of the effect of these parameters on the PHP performance are crucial to the optimum design of the PHP for various applications. Groll and his co-workers [2–4] considered a wide range of factors that affect the operating mechanism of PHP. Their work has vastly contributed to the knowledge of operational mechanism of PHP. Their review papers [5, 6] summarized the effects of the major parameters. To understand the flow characteristics of PHP, the flow visualization of PHP has been conducted by many researchers [7–9]. The early observation of oscillating and circulating motions were presented by Tong et al. [7]. Khandekar and Groll [8] observed the relationship between the magnitude of supplied heat and the pattern of motion. They reported that the flow motion changed from oscillating motion to circulating motion with the increasing of heat supply. The study of Xue and Qu [9] noted that the thermo-hydrodynamic characteristics and thermal performance have a direct relationship with the distribution of vapor bubbles and liquid slugs. Their investigations have brought a good progress in understanding PHP operational physics, although the cause of circulation and oscillation motions was yet unrevealed. Jang et al. [10] recently reported an experimental result of PHP in which they investigated the effect of nonuniform heating that they considered as the practical operating condition compared to a uniform heating. They charged HFE-7100 fluid in 1.8 mm-diameter tube PHP of eight turns. They observed different mode of flow direction for various heat inputs. Oscillation motion was dominant at low heat input and the motion changed to circulation pattern at higher heat input. As the heat input increased further, dryout occurred. Several theoretical and numerical studies have been performed to understand the heat transfer mechanism of PHP. Shafii et al. [11] developed one-dimensional modeling to predict thermal behavior and heat transfer of PHP. Following the Shafii’s model, many mathematical models have been proposed to analyze the complicated behavior of PHP motion [12–15]. In particula, Bae et al.

[15] investigated the effect of film dynamics on fluid motion and thermal performance of PHP. Their results showed a good prediction of thermal performance in comparison with experimental data, not only for a vertical PHP but also for horizontal and inclined PHP with various different parameters. In the computational fluid dynamics approach to the simulation of PHP, Lin et al. [16] presented a two-dimensional simulation based on volume of fluid (VOF) and mixture model. Their study showed that the mixture model is more suitable for the twophase flow simulation of PHP. Pouryoussefi and Zhang [17] also performed a two-dimensional simulation and reported the chaotic flow behavior in a multi-turn closed PHP. The formation of perfect vapor, liquid plugs, and the liquid film were observed in their study. Liu and Chen [18] simulated a flat-plate oscillating heat pipe (FP-OHP) using a three-dimensional model. They pointed out that the generation and characteristic of vapor slugs and liquid plugs come from the self-growth and coalescence of dispersed bubbles. Daimaru et al. [19] reported numerical simulation and experiment in an application of PHP to on-orbit systems. Although there have been many experimental and numerical work to investigate the operational mechanism and the thermal performance of PHPs, only a few numerical work have been carried out on predicting performance of PHP and yet successful simulation of actual pulsating actions seems very scarce. To develop a simulation model of PHP with better performance, it seems necessary to further investigate the flow patterns of working fluid in the evaporator and condenser sections. In the present study [20], a PHP was made of transparent Pyrex tubes for visualization of flow regimes. Both heating and cooling were provided by water jackets in which inlet and outlet temperatures were precisely measured and these were used to calculate heat transfer rates. The flow motion was recorded using a high-speed camera. A three-dimensional numerical simulation model has been developed using ANSYS Fluent. 2. Experiment 2.1. Experimental apparatus The test facility consists mainly of a heat pipe as the test section, data acquisition system, flowmeters, two circulating water baths. Fig. 1 shows a schematic diagram of experimental apparatus. Two circulating baths (JEIO TECH, RW-0525) were used to supply hot and cold water into heating and cooling water jackets. The hot and cold water were normally maintained at 80 °C and 25 °C, respectively. The water flow rate was variable but was normally 1 LPM and it was measured using a rotameter-type flowmeter. After the evacuating and charging process, the circulating baths flow the water through the chambers. The pulsating heat pipe was constructed using transparent Pyrex tubes of 3.02 mm of the outer diameter and 1.85 mm of inner diameter. Sixteen tubes were equally spaced by 26.2 mm apart and these make eight turns as shown in Fig. 2. The PHP was designed to operate in either closed or open type by opening or closing the connecting tube between the two ends. The major dimensions of the PHP are also found in Fig. 6 in the simulation part. The test section consists of four parts: adiabatic-1, adiabatic-2, heating, and cooling section as indicated in Fig. 6. Tests were conducted for vertical orientation with bottom heating. For working fluid, R123 was chosen based on the performance of various fluids reported in the prior work [3, 21]. The fluid motions were recorded using a high-speed camera (FASTCAM SA4). Twelve thermocouples were attached to the wall of adiabatic1 section as shown in Fig. 2(a). The inlet temperature of water jacket was measured with one thermocouple but the outlet temperature was measured with three thermocouples to get an aver-

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

3

Fig. 1. Schematic diagram of experimental apparatus.

Table 1 Experimental parameters.

Table 2 Slug velocity calculation.

Parameter

Value

From frame

To frame

From time

To time

t (s)

Velocity (m/s)

Tube inner diameter Number of tubes Orientation Water temperature: - Heating - Cooling Working fluid Filling ratio (FR) Water flow rate

1.85 mm 16 (8 turns) Vertical (bottom heating)

40 53 61 68 74

53 61 68 74 79

0.808 0.834 0.850 0.864 0.876

0.834 0.850 0.864 0.876 0.886

0.026 0.016 0.014 0.012 0.010

0.46 0.75 0.86 1.0 1.2

80 °C 25 o C R123 50–60% 1 LPM

aged value. Temperature data were monitored and recorded by a LabVIEW-based data acquisition system (National instruments, Ni SCX1 10 0 0). The thermal performance of the PHP was evaluated by measuring heat transfer rate, which was calculated based on the heat balance equation given below using the inlet and outlet temperatures of heating and cooling jackets and the water flow rate.

q = m˙ × C p × |(Tout − Tin )|

(1)

The system was designed with two 3-way valves installed on the adiabatic 2 section for vacuuming the PHP and charging the working fluid. Before charging R123, the PHP was evacuated to pressure less than 0.3 kPa by using the vacuum pump. Then, the working fluid was charged into the evacuated PHP with the preweighted amount of liquid for the desired filling ratio (FR). The experimental parameters are summarized in Table 1. 2.2. Experimental results A set of preliminary tests were conducted to observe the performance of the test PHP and to check the energy balance between the heating and cooling sections. Seven tests were chosen here to discuss on the flow pattern and heat transfer rate. These tests were divided into two groups based on filling ratio. Group A of 60% of filling ratio includes test no. 16, 17, 18 and 22 and group B of 50% of filling ratio includes test no. 19, 20 and 21.

2.2.1. Flow pattern visualization When a boiling process takes place, a variety of configurations known as flow patterns were observed. The particular flow pattern depends on the condition of pressure, flow velocity, heat flux, characteristics of fluid, and channel geometry. Fig. 3(a) and (b) shows the flow pattern of cooling and heating section. Flow patterns in the heating and cooling sections are generally different probably because of difference in heat transfer process of boiling and condensation. It can be seen that the major flow patterns in the heating section are annular, semi-annular, slug, and churn flow. In the cooling section, the major flow patterns are bubbly and annular flow. In the cooling section where the condensation occurs, the pressure difference between liquid slug and vapor plug is small. Therefore, the velocity of liquid slug and vapor plug is slow. The flow pattern representing this slow motion is slug slow. In the heating section, along the length of tube, the flow pattern is bubbly, slug, or annular flow. Flow visualizations were also used in estimating the velocity of fluid flow to determine whether the flow is laminar or turbulent. Fig. 4 shows the flow motion of slug flow in a bend of the heating section. The motion of slug flow from frame 40 to 79 of test 16 was analyzed. In this figure, the slug flow moves upward. The velocity of slug flow was calculated and presented in Table 2. The average velocity of fluid flow fluctuates in the range from 0.8 to 2 m/s. The Reynolds number calculated for liquid at this velocity range lies in the range from 6600 to 16,0 0 0, implying that the flow is turbulent. 2.2.2. Heat transfer rate The average heat transfer rate was calculated using Eq. (1) using the temperature data when the system reached steady state in

4

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317 Table 3 Summary of test results. Group

Test no.

Filling ratio (%)

Heat transfer rate (W)

Flow regime

State

A

16 17 18 22 19 20 21

60

180.7 189.6 169.0 173.3 181.6 186.1 171.4

Oscillating flow Circulating flow Circulating flow Oscillating flow Circulating flow Circulating flow Circulating flow

Stable Stable Reverse observed Stable Stable Stable Stable

B

50

Table 4 Measurement uncertainty of heat transfer rate. Variable

Uncertainty

Mass flow rate Temperature difference Heat transfer rate

± 2% ± 5% ± 5.4%

about one hour. The test results are summarized in Table 3. The oscillation motion was observed in test 16 and 22. The rest of the tests showed circular motion. The difference in heat transfer rate between two motions is small as seen in Table 3. The average heat transfer rate of group A is 181.2 W except test 18 because of reversal of motion which may cause the reduction of heat transfer. The average heat transfer rate of group B is 179.7 W. The uncertainty in the measurement of heat transfer rate was ± 5.4% as given in Table 4. The difference on the characteristics of fluid motion could be simply described with the direction of flow: in the case of circulating motion, the fluid normally travels in one direction, and in the case of oscillating motion, the fluid oscillates around a certain length. Circulating and oscillating motions could be identified by observing the wall temperature of two adjacent tubes. The main parameters resulting in different motion types are the filling ratio and the initial distribution pattern of the liquid and vapor throughout the heat pipe. The filling ratio is controlled in the

experiment and is given in Table 2. But the initial distribution pattern of liquid and vapor is random and it changes also as the test is performed repeatedly even with the same charged fluid. The circulating motion is normally observed when the initial distribution of the slugs and plugs is random and scattered in small size inside the pipe. Meanwhile, the oscillating motion normally occurs when the liquid is distributed as long slug columns. Such an effect of initial fluid distribution on flow behavior has been also reported by Tong et al. [7] and Daimaru et al. [19].

2.2.3. Wall temperature Fig. 5 illustrates the variation of wall temperatures in test 16, 18 and 19 for one hour after the steady-state was reached. Tests 16 and 19 represent for oscillating and circulating motions, respectively. The line numbers correspond to the location of twelve thermocouples as shown in Fig. 6. The odd-numbered and the evennumbered points are indicated by blue and red color respectively to express the characteristic of motion. As seen in the figure, in the case of oscillating motion, the temperature fluctuate around the average temperature of the heating and cooling water temperatures (about 56 ˚C). However, in the case of circulating motion, the odd-numbered and the even-numbered temperatures fluctuate at different level. Typically, in test 19, the odd-numbered temperature fluctuates around 51 ˚C and the even-numbered temperature fluctuates around 60 ˚C. The reason for this difference between two motions is the characteristic of fluid motion.

Fig. 3. Flow pattern of fluid flow at (a) heating section and (b) cooling section

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

In the circulating case, the fluid travels in one direction. Thus, the temperature at downstream points of the heating section is higher and the temperature at downstream points of the cooling section is lower. Otherwise, in the case of oscillation motion, the fluid moves back and forth over a certain length. Therefore, the temperature fluctuates in the range of average temperature. It also implies that the direction of motion can be determined by observing the temperature at the odd and even-numbered points. For example, in test 19, the temperature at the even-numbered points were higher than that of odd-numbered points. This indicates that the circulation direction of test 19 is counterclockwise. In test 18 of Fig. 5, the reversed direction of motion was observed repeatedly. This flow direction reversing behavior may degrade the thermal performance of PHP, but the reduction in heat transfer rate is small, about 10% as seen in Table 3. The reason for the reversing behavior seems to be non-uniform distribution of working fluid inside tubes. 3. Numerical simulation As the pulsating heat pipes have recently drawn a great attention in the field of heat dissipation technique for high and localized heat-source devices, there have been many theoretical studies in pursuit of computational simulation of pulsating actions of a PHP. The literature indicates, however, that successful simulation of actual pulsating actions seems scarce. In the present study, a threedimensional computational fluid dynamics modeling has been developed for pulsating heat pipe and the simulation has been able to show actual pulsating and circulating actions. The model the PHP consists of four sections: evaporator, condenser, adiabatic 1 and adiabatic 2 as show in Fig. 6. According to the previous experimental study [18], the filling ratio in the range of 50-60% for R123 is suitable to achieve high thermal performance. The present numerical simulations were performed for these two filling ratios: 50% and 60%. The cases FR = 50% and FR = 60% are named Case 1 and Case 2. Various parameters associated with the PHP geometry are found in Table 1. The present simulations used ANSYS Fluent v14.5.

5

low but enough to start the pulsating action, the dominant flow pattern of the working fluid is the slug flow. Under this condition, the oscillation of the working fluid is quite random with very low velocity. As the heat flux becomes higher, the flow pattern is likely to change to semi-annular flow and even the annular flow. From the experimental observation of flow patterns for the practical thermal conditions of 80°C wall temperature, the slug flow is the dominant flow pattern. Base on this conclusion, the Volume Of Fluid (VOF) two-phase flow model is a compatible method to simulate the PHP. The VOF model is a surface-tracking technique applied to a fixed Eulerian mesh. In the VOF method, the positions of vapor, liquid, and interface in the computational cells are represented by the volume fraction α v and α l , where subscripts v and l represent vapor and liquid, respectively. The liquid phase only exists in the cell where α v = 0 and the vapor phase only exists where α v = 1. Naturally, the vapor–liquid interface locates in the cell where 0 < α v < 1. In each control volume, the volume fractions of all phases sum up to unity,

αv + αl = 1

(2)

In general, we can specify the primary and secondary phases whichever way we prefer. In this model, vapor was defined as the primary phase to improve the solution stability and is treated as a compressible ideal gas. The tracking of the interface(s) between the phases is generally accomplished by the solution of a continuity equation for the volume fraction of one (or more) of the phases. For the phase q, this equation has the following form



1

ρq

 n  ∂ (α ρ ) + ∇ (αq ρqvq ) = Sαq + (m˙ pq − m˙ qp ) ∂t q q p=1

(3)

where m˙ qp is the mass transfer from phase q to phase p and m˙ qp is the mass transfer from p phase to phase q. By default, Sαq is the source term on the right-hand side of Eq. (3) and equal to zero. Momentum equation is as follows.

3.1. Governing equations

   ∂ (ρv ) + ∇ .(ρvv ) = −∇ p + ∇ . μ ∇v + ∇vT + ρ g + F ∂t

The previous study on visualization of PHP reveals that the flow pattern of the working fluid inside the PHP is generally related to the magnitude of applied heat flux [9]. When the heat flux is very

As discussed in the experimental results, the flow inside the PHP tube is likely to be turbulent, thus turbulent modeling is included in this study. The transport equations for k and ε in the

Fig. 4. Flow motion of slug flow in a bend of heating section.

(4)

6

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

Fig. 2. (a) Schematic of PHP, (b) photograph.

realizable k−ε model are as follow.

 ∂ ∂  ∂

μ  ∂k + Gk + Gb − ρε − ρ ku j = μ+ t (ρ k ) + ∂t ∂xj ∂xj σk ∂ x j YM + Sk

(5)

Fig. 5. Variation of wall temperatures of test 16, 18 and 19.   ∂ ∂  ∂

μ ρC2 ε 2 ε μ + t ∂ ε /∂ x j + ρC1 Sε − ρε u j = √ + C1ε C3ε Gb + Sε (ρε ) + ∂t ∂xj ∂xj σε k k + vε

(6)

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

7

by the presence of the component phases in each control volume,

ρ = αv ρv + (1 − αv )ρl

(11)

μ = αv μv + (1 − αv )μl

(12)

The ANSYS Fluent supports two models for simulating the interphase mass transfer through evaporation-condensation [22]. With the VOF formulation, only Lee model could be used. Lee [22] model is a simplified saturation model for evaporation and condensation processes. The key premise of this model is that phase change is driven primarily by deviation of interfacial temperature from Tsat and the phase change rate is proportional to this deviation. Therefore, the phase change occurs while maintaining temperatures of the saturated phase and interface equal to Tsat . The model assumes that mass is transferred at constant pressure and quasi-thermoequilibrium state according to the following relations. If Tl > Tsat (evaporation):

m˙ lv = coeff∗αl ρl

(Tl − Tsat ) Tsat

(13)

If Tv < Tsat (condensation):

m˙ vl = coeff∗αv ρv

Fig. 6. Configuration of the PHP numerical model.

where

C1 = max 0.43,

η 

η+5

1 2

,

η=S

k

ε

, S=



2Si j Si j

(7)

ρ k1 ∇ α2 ( ρ1 + ρ2 )

  ∂ (ρ E ) + ∇ .(v(ρ E + p) ) = ∇ . ke f f ∇ T + Sh ∂t

(9)

The VOF model treats energy E and temperature T as massaveraged variables.

αq ρq Eq q=1 αq ρq

q=1

E = n

(14)





q = h f Tw − T f + qrad

(15)

3.2. Geometry and mesh The computational domain was modeled to simulate the experiment given in the previous section, as shown in Fig. 6. The tube wall thickness was not modeled in the present analysis to avoid treating a conjugate heat transfer problem. To compare the tube wall temperature behavior between the simulation and the experiment, twelve points were chosen on the adiabatic region 1 as marked with a square symbol in Fig. 6. The meshing configuration of the domain is shown in Fig. 7. The tube lengthwise mesh size was 1.85 mm and the cross sectional mesh layout is also shown in the figure. The thickness of the first layer from the wall was 0.096 mm. The hexahedral mesh was used for entire domain and the total number of elements was 150,570.

(8)

where ρ is the volume-averaged density computed and the above equation shows that the surface tension source term for a cell is proportional to the average density in the cell. Energy equation is as follows.

n

Tsat

When a constant temperature boundary condition is applied to the wall, the heat flux to the wall from a fluid cell is computed as

In these equations, Gk represents the generation of turbulence kinetic energy due to the mean velocity gradients. Gb is the generation of turbulence kinetic energy due to buoyancy. YM represents the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate. C2 and C1ε are constants. σ k and σ ε are the turbulent Prandtl number for k and ε, respectively. Sk and Sε are user-defined source terms. The current model used the Fluent default settings for these parameters. The details of the model is found in ANSYS Fluent documents [22]. The continuum surface force (CSF) model in ANSYS Fluent [22] has been implemented such that the addition of surface tension to the VOF calculation results in a source term in the momentum equation. The general form of surface tension equation is the following.

Fvol = σ12

(Tsat − Tv )

(10)

where Eq for each phase is based on the specific heat of that phase and the shared temperature. The properties ρ and keff (effective thermal conductivity) are shared by the phases. The source term, Sh , contains contributions from radiation, as well as any other volumetric heat sources. The fluid properties in the above governing equations are determined

3.3. Initial and boundary conditions The initial and boundary conditions were set as close to experimental conditions as possible to simulate the experiment given in the previous section, as was the domain geometry. A random liquid-vapor distribution was initialized as shown in Fig. 8. The initial temperature of both liquid and vapor were set to 52.45 °C and the initial pressure was set to the saturation pressure of initial temperature. The constant temperature boundary conditions were applied to the wall boundary of condenser and evaporator. The condenser wall was set to 25 °C and the evaporator wall was set to 80 °C. The saturation temperature of working fluid was provided as a function of pressure. A set of preliminary simulations revealed that weather the density was set to constant or to variable was important. Thus vapor density was calculated by the ideal gas law and liquid density was given by a function of temperature. The major parameters of the simulation including the boundary conditions are summarized in Table 5. The equations of saturation temperature and liquid density of R123 were obtained by least-square

8

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

Fig. 7. Illustration of meshing on tube cross section and lengthwise.

Table 5 Parameters and boundary conditions for numerical simulations. Parameters

Value

Working fluid Tube inner diameter Wall boundary conditions Heating section Cooling section Filling ratio Initial conditions Pressure Temperature

R123 1.85 mm

3.4. Results and discussions

80 °C 25 °C 50%, 60% 212.46 kPa 323.15 K

fitting of the properties from the REFPROP 9.0 and these are as follow.

In the present study, a three-dimensional computational fluid dynamics modeling has been developed for a closed-type pulsating heat pipe and the simulation has been able to show actual operation of the PHP. In order to validate the proposed numerical model, the simulation results were compared with the experimental measurements. In the present simulations, only the circulating motion was observed in both cases of different filling ratio. The volume fraction distribution and flow direction, wall temperature, and heat transfer rate are discussed. The total fluid mass inside the heat pipe was checked in every calculation by exporting the integration of all phase mass in Fluent and it was found that the total mass was conserved throughout the transient calculation. The initial liquid-phase fluid filling in the heat pipe was typically 50% or 60% so there was free vapor volume to accommodate the change of density due to variable density modeling.

Tsat =268.269+3.6799×10−4 P − 6.1422×10−10 P 2 +4.3892×10−16 P 3 (16)

ρL =978.327+13.06T −0.0759T 2 +1.70539×10−4 T 3 −1.493×10−7 T 4 (17)

3.4.1. Volume fractions and circular motion Fig. 9 illustrates volume fractions of liquid and vapor at different times (red color represents the vapor and blue color represents the liquid). To show the transient characteristic of the motion, the void fraction contours are given together for six different times of 0, 1.5, 3.4, 5.6, 6.3 and 25 s. It is observed that long liquid columns are produced by condensation, broken into smaller col-

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

Fig. 8. Random liquid–vapor distribution at the initial time (t = 0) for the case of FR = 50%. Table 6 Fluid velocity calculation of simulation results From time (s)

To time (s)

Distance (mm)

t (s)

Velocity (m/s)

3.56 3.58 3.60 3.62 3.64

3.58 3.60 3.62 3.64 3.66

19 18 17 16 16

0.02 0.02 0.02 0.02 0.02

0.95 0.9 0.85 0.8 0.8

umn by evaporation, and they are circulating along the tubes, in most cases, in the counterclockwise direction. At the initial time t = 0 s, the vapor plug and liquid slug distribute randomly inside the tube. At time t = 1.5 s, the start-up process of PHP begins. The liquid slugs combine to a long liquid column locating on the cooling section. At time t = 3.4 s, because of the difference in pressure between the heating and cooling regions inside the tube, the liquid columns start to travel in one direction. As we can see, in this case, the liquid columns move in the counterclockwise. At time t = 5.6 s, the liquid columns still move downward. However, the liquid columns move slower and start to reverse the direction of motion. At time t = 6.3 s, the liquid column motion is reversed to the clockwise direction. These liquid columns combine with the smaller liquid slugs ahead. This behavior of liquid column motion occurs repeatedly. The fluid velocity in the simulation was sampled and given in Table 6 to compare it with the experimental data shown in Table 2. The fluid velocity in the experiment ranged from 0.46 m/s to 1.20 m/s and in the simulation it ranged from 0.8 m/s to 0.95 m/s, indicating reasonable agreement in fluid velocity. The motion pattern in the experiment was mainly circulating, but sometimes oscillating pattern was also observed. But in the simulation, the motion pattern was circulating in all cases. 3.4.2. Wall temperature In the experimental measurement, the wall temperature is normally used to evaluate the motion behavior of working fluid. It

9

represents the temperature of fluid inside the tube at the time the fluid passes the measuring position and this wall temperature could be more strongly affected by the liquid than the vapor because the thermal conductivity of liquid is much higher than that of vapor. Fig. 10 shows the predicted wall temperatures for the case of 60% filling ratio and those of 50% filling ratio also look similar to these. The locations of twelve points of temperature measurement are indicated in Fig. 6. The fluctuating pattern of the calculated wall temperatures showed good agreement with the experimental observation shown in Fig. 5. The wall temperature variations are characterized in two distinct patterns; one for odd-numbered points (blue lines) and the other for even-numbered points (red lines). The wall temperatures at odd-numbered points fluctuate in the range which is lower than the range of even-numbered points. In general, oscillating pattern shows smaller temperature difference but circulating pattern brings large temperature difference. So the pattern in Fig. 10 can be compared with Test 19 in Fig. 5, both circulating patterns. This unique characteristic of wall temperature behavior can be attributed to the characteristic of the circulation motion. The temperature of the fluid which leaves the heating section is generally higher than the temperature of the fluid which leaves the cooling section. Since the direction of circulation in this case is counterclockwise, the temperature at the even-numbered points are always higher than that of the odd-numbered points. This feature can imply that by observing the wall temperatures, the regime of motion and the direction of motion can be determined. The figure also shows the difference in the amplitude of wall temperature fluctuation between the odd-numbered and the even-numbered points. Such a difference in fluctuation amplitude was also observed in the experiment. 3.4.3. Heat transfer rate The heat transfer rate from the wall to the fluid can be obtained by integrating the heat flux across the wall. The integration was done separately for the heating and the cooling sections. Fig. 11 illustrates the time-variations of heat transfer rates in both heating and cooling sections for the two different filling-ratio cases. In the early period, the heat transfer rates of heating and cooling fluctuate with large amplitude. It looks that the fluctuation of heating rate is more severe than the cooling rate. However, this large fluctuation soon disappears after about ten seconds and the heating rate becomes equal to the cooling rate, indicating quasi-steady state operation of the PHP. The average heat transfer rates were calculated by integrating the value from t = 20 s to t = 25 s in which quasi-steady state operation is achieved. The results are compared with the experimental data in Fig. 12. Although the agreement is fairly good, it is noted that the experimental data showed that the average heat transfer rate was about 180 W and there was little difference between the two filling ratios. Meanwhile, the simulation results show that the heat transfer rate of 60% filling ratio is about 10% higher than that of 50% filling ratio. The simulation shows higher heat transfer rate in FR = 60% than that of FR = 50%, while the experimental results show little difference between the two filling ratios. Theoretically, higher filling ratio provides more liquid to evaporate thus results in higher heat transfer rate. But too high filling ratio provides little room to accommodate the vapor to condense in cooling section. Therefore it can be said that the simulation shows more correct trend of the effect of filling ratio on heat transfer than the experimental data, but the difference seems small that it is hard to neglect any counterbalancing mechanisms in heat transfer. It is generally known that the charge ratio ranging from 20% to 80% will make the PHP device operate as a true pulsating heat pipe. An optimal

10

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

Fig 9. Transient behavior of volume fractions of liquid (blue) and vapor (red) (FR = 60%).

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

Fig. 10. Variation of wall temperatures (FR = 60%).

11

Fig. 12. Comparison of heat transfer rate between experiment and simulation.

the first cell, the y-plus of 30 was chosen. From the experiment, it was observed that the magnitude of velocity was in the range of near 1.2 m/s. With this velocity and the y-plus of 30, the calculated layer thickness was 0.0906 mm and the actual mesh value was 0.0961 mm. One case with smaller first-layer of 0.059 mm was tested and the predicted heat transfer rate increased by 13% and thus it overpredicted the experimental data by 18%.

4. Conclusion

Fig. 11. Heat transfer rates of heating and cooling: (a) FR = 50% and (b) 60%.

charge ratio exists for each particular PHP setup. In our study, the optimum charge ratio is around 50–60%. The grid independence was checked with the focus on the wall function point of view. In this study, the standard wall function was applied for the realizable k−ε turbulence model. For the standard wall function, the y-plus value should be typically in the range of 30 and 300. To determine the distance from the wall for

An experimental and analytical study of pulsating heat pipe (PHP) has been carried out in order to develop a more reliable numerical simulation model. The test PHP was made of transparent Pyrex tubes for visualization purpose. The inner diameter of tube was 1.85 mm and a total of sixteen tubes formed eight turns of PHP. Both heating and cooling were provided by water jackets in which inlet and outlet temperatures were precisely measured and these were used to calculate heat transfer rates. The working fluid was R123. Visualization using a high-speed camera showed various flow patterns as well as the fluid motions. The wall temperature measured at various locations revealed its relation to the fluid motion and the direction of circulation. The circulation motion was dominant in most tests. The heat transfer rate was measured and the difference for filling ratio of 50 and 60% was little. A three-dimensional computational fluid dynamics modeling has been developed for pulsating heat pipe using ANSYS Fluent v14.5. The VOF model with variable density and vapor pressure relation successfully simulated the circulating motions of PHP. The simulation of the experiment conducted in the present study showed the same circulating motion as the experiment. The predicted wall temperatures also showed the same indication of mode of flow motion and the flow direction as observed in the experiment. The predicted heat transfer rates agreed well with the experimental data within 5%. The success of the simulation of the experiment implies that the realizable k-ε turbulence model is appropriate to PHP simulation and the use of variable density for liquid and vapor and the vapor pressure equation is crucial.

Declaration of Competing Interest I have no conflict of interest in publishing the manuscript titled “An experiment and three-dimensional numerical simulation of pulsating heat pipes” in the International Journal of Heat and Mass Transfer.

12

D.-T. Vo, H.-T. Kim and J. Ko et al. / International Journal of Heat and Mass Transfer 150 (2020) 119317

CRediT authorship contribution statement Duy-Tan Vo: Investigation, Writing - original draft. Hyoung-Tak Kim: Visualization, Validation. Junghyuk Ko: Methodology, Investigation. Kwang-Hyun Bang: Conceptualization, Supervision, Writing - review & editing. Acknowledgment This work has been supported by the National Research Foundation of Korea (Grant No. 2018R1D1A1B07049027). Supplementary materials Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.ijheatmasstransfer. 2020.119317. References [1] H. Akachi, Structure of a heat pipe, U.S. Patent No. 4921041 (1990). [2] P. Charoensawan, S. Khandekar, M. Groll, P. Terdtoon, Closed loop pulsating heat pipe Part A: parametric experimental investigations, Appl. Therm. Eng. 23 (20 03) 20 09–2020. [3] S. Khandekar, N. Dollinger, M. Groll, Understanding operational regimes of closed loop pulsating heat pipes: an experimental study, Appl. Therm. Eng. 23 (2003) 707–719. [4] H. Yang, S. Khandekar, M. Groll, Operational limit of closed loop pulsating heat pipes, Appl. Therm. Eng. 28 (2008) 49–59. [5] Y Zhang, A Faghri, Advances and unsolved issues in pulsating heat pipes, Heat Transf. Eng. 29 (2008) 20–44. [6] X. Han, X. Wang, H. Zheng, X. Xu, G. Chen, Review of the development of pulsating heat pipe for heat dissipation, Renew. Sustain. Energy Rev. 59 (2016) 692–709.

[7] B.Y. Tong, T.N. Wong, K.T. Ooi, Closed-loop pulsating heat pipe, Appl. Therm. Eng. 21 (2001) 1845–1962. [8] S. Khandekar, M. Groll, An insight into thermal-hydrodynamic coupling in closed loop pulsating heat pipes, Int. J. Therm. Sci. 43 (2004) 13–20. [9] Z.H. Xue, W. Qu, Experimental and theoretical research on an ammonia pulsating heat pipe: new full visualization of flow pattern and operating mechanism study, Int. J. Heat Mass Transf. 106 (2017) 149–166. [10] D.S. Jang, H.J. Chung, Y. Jeon, Y. Kim, Thermal performance characteristics of a pulsating heat pipe at various nonuniform heating conditions, Int. J. Heat Mass Transf. 126 (2018) 855–863. [11] M.B. Shafii, A. Faghri, T. Zhang, Thermal modeling of unlooped and looped pulsating heat pipes, J. Heat Transf. 123 (2001) 1159–1172. [12] R. Senjaya, T. Inoue, Oscillating heat pipe simulation considering bubble generation Part I: Presentation of the model and effects of a bubble generation, Int. J. Heat Mass Transf. 60 (2013) 816–824. [13] H. Peng, P.F. Pai, H. Ma, Nonlinear thermomechanical finite-element modeling, analysis and characterization of multi-turn oscillating heat pipes, Int. J. Heat Mass Transf. 69 (2014) 424–437. [14] I. Nekrashevych, V.S. Nikolayev, Effect of tube heat conduction on the pulsating heat pipe start-up, Appl. Therm. Eng. 126 (2017) 1077–1082. [15] J. Bae, S.Y. Lee, S.J. Kim, Numerical investigation of effect of film dynamics on fluid motion and thermal performance in pulsating heat pipes, Energy Convers. Manag. 151 (2017) 296–310. [16] Z. Lin, S. Wang, R. Shirakashi, L.W. Zhang, Simulation of a miniature oscillating heat pipe in bottom heating mode using CFD with unsteady modeling, Int. J. Heat Mass Transf. 57 (2013) 642–656. [17] S.M. Pouryoussefi, Y. Zhang, Analysis of chaotic flow in a 2D multi-turn closed-loop pulsating heat pipe, Appl. Therm. Eng. 126 (2017) 1069–1076. [18] X. Liu, Y. Chen, Fluid flow and heat transfer in flat-plate oscillating heat pipe, Energy Build. 75 (2014) 29–42. [19] T. Daimaru, H. Nagai, M. Ando, K. Tanaka, A. Okamoto, H. Sugita, Comparison between numerical simulation and on-orbit experiment of oscillating heat pipes, Int. J. Heat Mass Transf. 109 (2017) 791–806. [20] V.D. Tan, An Experimental and Analytical Study on Pulsating Heat Pipe, Korea Maritime and Ocean University, 2019 M.S. Thesis. [21] H.-T. Kim, H.-K. Park, K.H. Bang, Heat dissipation of sealed LED light fixtures using pulsating heat pipe, J. Korean Soc. Mar. Eng. 36 (2012) 64–71. [22] ANSYS Fluent Theory Guide, v.14.5 (2011).