- Email: [email protected]
Numerical and experimental investigations of electronic evaporative cooling performance with a coiled channel
Applied Thermal Engineering 94 (2016) 256–265
Contents lists available at ScienceDirect
Applied Thermal Engineering j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / a p t h e r m e n g
Research Paper
Numerical and experimental investigations of electronic evaporative cooling performance with a coiled channel Xiang Yin a, Feng Cao a,*, Lei Jin a, Bin Hu a, Pengcheng Shu a, Xiaolin Wang b a b
School of Energy and Power Engineering, Xi’an Jiaotong University, 28 Xianning West Road, Xi’an, 710049, China School of Engineering and ICT, University of Tasmania, Private bag 65, Hobart, TAS 7001, Australia
H I G H L I G H T S
• • • • •
Two-phase flow cooling performance of a moving electronic unit was studied. The distributions of the vapor and the liquid phase in a 3D snake pipe were simulated. An experimental setup was established to verify the simulation results. The heat concentration in the curve of the channel affects the cooling effect directly. An optimized method was developed to improve the overall cooling performance.
A R T I C L E
I N F O
Article history: Received 12 July 2015 Accepted 22 October 2015 Available online 3 November 2015 Keywords: Electronic evaporative cooling Numerical simulation Coiled channel VOF Phase distribution Cooling performance
A B S T R A C T
Numerical and experimental investigations on two-phase evaporative cooling performance for a moving electronic unit with high heat flux have been performed in a coiled channel heat exchanger. Based on the VOF model, the phase distribution of a 3D coiled channel was investigated. It was found that the main vapor streamline flowed through the inner side of one bend channel to the next one, which increased the contact frequency of the vapor and the liquid, improving the cooling performance. Moreover, the results were also verified by the experiment. The cooling effect was directly affected by heat concentration in the bend channel. An optimized method was implemented for this phenomenon. Decreasing the effective flow area and increasing the local velocity could establish a new force balance for the vapor and liquid, eliminating the heat concentration and improving the overall cooling performance. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction With the development of human civilization, electronic equipment sets, such as computers, radios, data servers and electric vehicles, intrude all respects of our daily lives. Some particular electronic equipment sets with smaller size have to deal with a large heat flux. Most of them need to work in relatively stable conditions with a small temperature fluctuation range. When the temperature exceeds certain values, the performance of the components goes down sharply, eventually stopping the device from working. Therefore, the cooling performance of electronic components is always an interesting topic for researchers. However, traditional cooling systems can no longer satisfy the requirements for the cooling performance of electronic components because of their limited cooling capability. Nowadays, two-phase evaporative cooling systems are widely applied in the electronic cooling field, e.g. cooling of radio [1], computers [2,3], or chips [4,5]. To improve
* Corresponding author. Tel.: +86 29 82663583; fax: +86 29 82663583. E-mail address: [email protected] (F. Cao). http://dx.doi.org/10.1016/j.applthermaleng.2015.10.127 1359-4311/© 2015 Elsevier Ltd. All rights reserved.
its application, the two-phase boiling performance should be further studied via experiments or numerical simulations. Extensive experimental literature reports on evaporative cooling are available, while only a few deal with evaporative cooling in the macroscopical evaporator by using numerical methods. In a boiling flow, both the convective boiling and the nucleate boiling regimes affect heat transfer. Different mechanisms constitute the main regime in different flows with different shapes [6]. Aixiang Ma et al. [7] concluded that the inner geometry of evaporators has little influences on the two-phase boiling heat transfer when heat loads are high [8]. Moreover, the Reynolds number has a direct effect on heat transfer coefficient of forced convective boiling [9]. In addition, vapor friction has a little influence on the heat transfer coefficient, while the heat flux and pressure have a stronger influence on the heat transfer [10]. That is to say, the effect of inner geometry on the boiling heat transfer is determined by the application and channel structure. A. Feldman et al. [11] studied the mechanism of boiling heat transfer and figured out that the heat transfer coefficient is independent of mass velocity and quality in the convective regime, but is affected by mass velocity and quality in the nucleate regime. In summary, the above mentioned papers
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
mainly showed that the two-phase cooling system could be better employed in the cooling of electronic components. Nevertheless, much research is still needed in view of better applications, especially in larger electronic units. Experimental studies on boiling flow are sometimes limited by practical conditions. It is hard to optimize the structure of the evaporator channel, if the boiling performance cannot meet the requirements. Thus, numerical simulations can be an alternative way to solve the problems. Some efforts in the development of twophase simulations, using the VOF model, have been done to trace the phase interface in a boiling flow. Woorim Lee et al. [12,13] successfully employed the VOF model to analyze the effect of cavity diameter and surface modification on boiling heat transfer. Moreover, numerical simulations have been widely used to investigate the heat transfer performance on the microstructure in a boiling flow [14,15]. The diameter and the number of cavities, as well as the structure of channels, play different roles in different boiling flows [16,17]. The numerical method could be well used to develop the surface configuration for boiling flows. At larger length scales, Christian Kunkelmann and Peter Stephan [18] conducted a numerical investigation on boiling of HFE7100 according to the experiment of Wanger et al [19]. They found that the departure diameter of a bubble is 2.3 mm, and the results are around 20% above the values obtained experimentally. However, theirs was a 2-dimensional simulation, which could not be directly applied into practice. It is worth pointing out that Y.W. Kuang et al. [20] simulated the flow boiling behavior in 3-dimensional pipes, whose diameters were 32 mm and 65 mm, respectively. The model could well simulate the flow pattern evolution, and capture the hydrodynamic as well as thermal mechanisms. But the geometric model is just a straight pipe. In addition, Z. Yang et al. [21] simulated the boiling flow in a coiled tube, and the results were in good agreement with the phenomenon observed in experiments. It is reasonable to claim that in those works the VOF model could be widely used to analyze the phase distribution. However, they did not analyze the correlation between the heat transfer performance and the phase distribution. As a whole, in the field of boiling flow, a number of investigations have been conducted with the VOF model, including the behavior of the bubble, the evolution of the flow pattern, and the performance of heat transfer. However, most of them are limited to a 2-d scale or to an ideal model. Actually, the influence of the bend on flow pattern or heat transfer of two-phase evaporative cooling has been neglected in the past [22]. The numerical investigation of the bend will be highlighted in this paper. In the present paper, an attempt was made to investigate the cooling performance difference of a coiled channel and the vapor behavior, which may affect the cooling performance. A numerical simulation, using the VOF model, and the corresponding experiments were conducted to study the evaporative cooling performance for electronic components in a horizontal coiled heat exchanger with a rectangular channel, which could increase the effective heat transfer area. The heat flux of every component ranges from 60 to 120 W/ cm2. The phase distribution was also studied. Numerical analyses of the cooling performance were verified by experimental data directly and indirectly. Moreover, the channel structure was optimized. 2. Simulation method 2.1. Model of simulation The tracking of phase interface, in the VOF model, was accomplished by solving the continuity equations for the volume fractions of different phases. When the number of the cells is large enough, the errors due to the difference between curves and straight lines could be neglected. The straight line could replace the smooth interface of the two-phase (see Fig. 1). To avoid truncation errors, which would cause a non-conservation of mass velocity between inlet and
257
Fig. 1. Interpolation of interface in VOF.
outlet, the mesh step-size must be selected to match the vapor size. However, in a transient simulation, the time-step is much sensitive to the mesh step-size and must be decreased to confirm stability and convergence, which usually cost much more time. In particular, in 3-dimensional simulations, the simulation period would be too long. Therefore, the mesh step-size should be increased at the expense of simulation accuracy, if the non-conservation is less than 5%. That is to say, some small-sized vapors are neglected.
2.2. Governing equations The fluid in the control volume is governed by three conservation laws: conservation of mass, conservation of momentum and conservation of energy. In the VOF method, the summation of the volume fraction of liquid and vapor equated to 1. The conservation of mass is:
∂α l S + ∇ (v l ⋅ α l ) = m ∂t ρl
(1)
∂α g S + ∇ (v g ⋅ α g ) = m ∂t ρg
(2)
where Sm is the mass source term that reflects the mass transfer between the two phases in a boiling flow, and l, v represent liquid and vapor, respectively. The conservation of momentum is:
∂ ( ρv ) + ∇ ⋅ ( ρvv ) = −∇p + ∇ ⋅ [ μ (∇v + ∇v T )] + ρg + Ftotal ∂t
(3)
The force, Ftotal, is the sum of all forces that affect the generation and movement of vapor, including the surface tension, the force between two phases and so on, excluding gravity. Here, the surface tension is highlighted because of phase changes. In a CSF model [23], the surface tension contributes to surface pressure, generating a force that affects the momentum of the fluid. The force follows that
Fg = σ k
∇ρ ρ ρl − ρ g 1 2 ( ρl + ρ g )
(4)
where k is the curvature of the interface, formed by the liquid and the vapor. It is defined as
k=
Δρ ∇ρ
(5)
258
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
The equation for conservation of energy is:
∂ ( ρE ) + ∇ ⋅ [v ( ρE + p )] = ∇ ⋅ (k eff ∇T ) + S E ∂t
(6)
where SE is the energy source term that satisfies the following equation:
S E = S V ∗ Δh
(7)
where Δh is the enthalpy difference between the liquid and vapor at the state point of saturation temperature in the condition of interface. 2.3. Mass and heat transfer theory
Fig. 2. Geometrical configuration.
In a boiling flow, when the liquid temperature increases above the saturation temperature or the vapor temperature decreases below the saturation temperature, phase change occurs, following mass and heat transfer through the interface of the liquid and vapor to obtain the thermodynamic equilibrium. Based on the kinetic theory, the Hertz Knudsen equation was proposed to calculate mass transfer through the free interface in 1953 [24].
Gf = β
M ( p * − p sat ) 2π RTsat
(8)
where β is the coefficient that reflects the ratio of absorbed vapor to the total amount of vapor entering the liquid surface. The flow was treated as in equilibrium at each time step. The vapor, which is generated from the heat wall, is considered as a constant size globular. Combined with the Clausius Clapeyron equation, the source term per volume (Sm) could be described by the flowing equation:
S m = G f ∗ A c V c = G f ∗ ( πd 2 ) ( 4 πr 3 3) = γ ∗ α g ρg
γ =
( T − Tsat )
⎛ ρl ⎞ 3 M β Δh ⎜ r 2π RTsat ⎝ ρl − ρg ⎟⎠
Tsat
(9)
(10)
3. Experimental facilities
The saturation temperature was treated as a constant number due to the small pressure drop less than 0.1 bar. In a control volume, when the temperature increases above the saturation temperature, some liquid transforms into vapor, to maintain a certain temperature. Conversely, some vapor transforms into liquid to keep the equilibrium. Thus, thermodynamic equilibrium is achieved. The control equation of phase change can be written as follows: If Tc > Tsat, the liquid transforms into vapor:
S m(l − g ) = γ ∗ α l ρl
(T − Tsat ) Tsat
to the fact that such working condition constitutes the worst case in the presence of a rotational movement that changes the position of the cooling equipment. The electronic components were replaced by thermal resistances with the same power that were attached to the double layer heat expanding device in order to decrease the heat flux and improve the cooling performance. Fig. 2 shows the configuration of the electronic components on the cooling heat exchanger with a rectangular coiled channel. The working fluid was R134a, and the secondary boundary was used on the electronic components. To improve the precision and reduce the simulation period, a structured grid was adopted, due to the absence of tortuous surfaces that would require unstructured grids to meet the bodyfitted performance. The mesh step was 0.5 mm in the channel and 0.25 mm near the boundary. The mesh step of the solid flied was 1 mm. In the case of transient simulation, using PISO for pressure interpolation allows lager relaxation factors, decreasing the simulation period. Time step was adapted to the global Courant number to obtain a stable simulation, using values of the order 1E-6 s. The k-ε model was adopted to compute turbulence effects. Additionally, the geometric reconstruction method was used to calculate the phase interface.
To demonstrate the reliability of numerical analysis, experimental investigations were conducted. The experimental equipment is shown in Fig. 3. The mass flow rate of the test heat exchangers was supplied by a sliding-vane working pump. The hand value could increase the adjustable range of the mass flow rate, besides
(11)
Else Tc < Tsat, the vapor transforms into liquid:
S m ( g −l ) = γ ∗ α g ρ g
(T − Tsat ) Tsat
(12)
where γ is calculated according to equation (10). In a certain simulation, γ is constant due to the fact that parameters like ρ, β, and d are definite value. When it comes to R134a, γ was set to 500 or more if the heat flux was too high to maintain thermodynamic equilibrium. 2.4. Geometrical configurations A horizontal working condition was performed to investigate the cooling performance for a moving electronic unit, according
Fig. 3. Experimental facility.
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
adjustment with inverter control. An electric heater, in front of the cooling heat exchanger, could control the vapor quality, which flows into the heat exchanger. The reservoir could improve the stability of the system. Four pressure transducers were placed at the inlet and outlet of the cooling heat exchanger and the pump. Meanwhile, twelve thermocouples were installed near the inlet and outlet of the heat exchanger and every electronic component. 4. Results and discussion 4.1. Phase distribution and cooling performance With the VOF model, the phase distributions could be observed directly, as shown in Fig. 4. The red color represents vapor, and the blue represents liquid. A steady phase distribution was selected to investigate the cooling performance. Although the flow pattern was not steady due to the changing flow parameters (e.g. pressure) in boiling flow, it had little effect on the overall flow tendency when the flow pattern developed well. Therefore, the cooling performance could be successfully investigated from a steady result. The
Fig. 4. Phase distribution of the whole channel.
259
phase distributions present here were selected at a relatively steady flow time. Vapor quality increased along the flow direction, and the amount of liquid decreased. The vapor mainly occupied the top side of the channel, as a consequence of gravity. Originally, the channel was full of liquid. When the electronic components were heated, the heat layer formed. Then, the vapor appeared where the temperature was higher than the saturation temperature. Fig. 4 shows that the vapor tended to gather at the inner side of the bend. Then, in a continuous bend, the main vapor streamline flowed from the inner side of one bend to the inner side of the next one, forming an angle between the main vapor streamline and the channel axis. Meanwhile, the continuous bend made the main vapor streamline change sharply, disturbing the phase interface, which should otherwise be smooth. Fig. 5a and b shows a 3-d translucent main vapor streamline map. The phase interface in the bend is rougher than the one in the straight channel. Thus, heat and mass occur through the whole phase interface form all directions, instead of a single surface. Furthermore, Fig. 5c and e shows the phase interface and the contact action of the liquid and vapor in the continuous bend and straight channel, respectively. Finally, Fig. 5d and f represents the section cut from planes S1 and S2. In the straight channel, the phase interface was smooth under low mass flow condition. Contact between liquid and vapor occurred just below the phase interface. However, in the continuous bend, the rough phase interface and the twisty main vapor streamline improved the contact opportunity of the liquid and the vapor. Mass and heat transfer occurred through the phase interface from all directions, thus increasing the contact frequency of the vapor and the liquid, and improving the cooling performance. Actually, the vapor converging tendency was not clearly observable, because new vapor was continuously generated at the top of the channel, slightly hiding the vapor tendency. Fig. 6a shows the phase contour upon getting rid of the top of the channel, and Fig. 6b represents the sketch map. It is easy to observe that the main vapor streamline does not flow at a constant section. The vapor always gathered at the inner side of the bend. If this vapor could flow downstream in time, the cooling performance would be enhanced; otherwise, a boiling crisis might occur. In our case, the vapor could flow successfully. The cooling performance improved at the continuous bend channel. However, the flow behavior was different at the last bend. The vapor flowed in a smoother manner.
Fig. 5. Difference flow performance in bend and straight channels.
260
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
Table 1 Working conditions of the experiment.
Group 1 Group 2
(a) Phase contour of the continuous
Mass flux/t/h
Total power of components/W
Air temperature/K
0.2~1.6 per o.2 0.2~1.6 per o.2
1000 1600
298 298
tions cut from the direction of SA, SB, and SC in Fig. 4. The main vapor streamline became thin or disappeared, while it was smooth in the channel after the last bend. The flow showed a different behavior, in agreement with the analysis above. Except the last bend, the vapor did not flow at a constant section. The vapor flowed at a constant section of the channel after the last bend, but the phase interface after the last bend was non-turbulent. The contact frequency between the vapor and the liquid was smaller. The cooling performance turned out to be worse than the others. Here, all the flow parameters were the same, except the vapor layer and the phase distribution. The vapor layer might cause a gradual decrease in the cooling performance, but different flow distributions would cause the cooling performance to decrease sharply. Consequently, the cooling performance of the last bend must be much worse than the others. 4.2. Experimental validation
(b) Sketch map of the vapor mainstream Fig. 6. Sketch map of the vapor flow through the continuous curve.
The contact frequency and area at the last bend are expected to be smaller. Therefore, the heat transfer and cooling performances would be worse at this bend than at the others in the same condition (i.e. same velocity, heat flux and heat area, geometrical surface). In addition, it was interesting to find that the vapor flow was not continuous at a constant section. Fig. 7b and c shows the sec-
(a) Phase contour of section SA
(b)
Phase contour of section SB
To demonstrate the analysis above, two groups of experiments were conducted (see Table 1). The electronic components from inlet to outlet were defined as numbers 1, 2, 3…, respectively, and numbers 5, 6, 8, and 10, which had similar heat flux, were selected to analyze the cooling performance. The cooling performance was characterized by the temperature of electronic components: the higher the temperature, the worse cooling performance, and vice versa. 4.2.1. Comparison of experiments and simulations Fig. 8a and b shows the electronic temperature discrepancy between experiment and simulation, whose total heat power is 1000 W and 1600 W, respectively. A similar tendency was found. The electronic component named number 1 might be characterized by a single phase flow, leading to a smaller discrepancy. When the phase changed, the error increased. The error values ranged from 1.23 K to 6.15 K, the largest one taking place under number 3. This may have happened because the subcooled temperature of the simulation is not the same as the experiment, leading to a smaller vapor quality under number 3 than the experiment. The reason why the cooling performance of the simulation was better than the experiment is that the actual thermal contact resistance between the electronic component and the cooling heat exchange was not as ideal as the simulation. The constant globular hypothesis and the simplifications of the VOF model, such as the repetitive computation of the interface advancement, also increased the error. 4.2.2. Indirect validation of the flow pattern The actual flow pattern in the channel was difficult to discover. Here, the flow performance could be deduced from the analysis of the cooling performance. According to the previous analysis, the cooling performance of number 10 was expected to be the worst. The convection heat transfer equation can be written as
Φ = hA ( Tw − Tf
)
(13)
while the heat conduction equation can be written as
(c) Phase contour of section SC Fig. 7. Phase contour of different section.
Φ ( T − Tw ) =λ e δ A
(14)
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
338 336
experiment
0.3
simulation
0.25
ΔNh6
0.2
ΔNh8
0.15
ΔNh10
ΔNh
temperature/K
334 332
0.1
330
0.05
328
0
326 0
2
4
number
6
8
364
(a)
ΔNh6
0.2
Simulation
360
Comparsion of cooling performance under 1000W
0.25
Experiment
362
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 mass flux /t*h-1
10
(a) Temperature comparison under 1000W
ΔNh8 ΔNh10
0.15
358
ΔNh/W*K-1
temperature/K
261
356
0.1
354
0.05
352 350 0
2
4 6 number
8
10
0 0.2
0.4
0.6 0.8 1 1.2 -1 mass flux /t*h
1.4
1.6
(b) Temperature comparison under 1600W (b) Comparison of cooling performance under 1600W
Fig. 8. Temperature comparison of experiment and simulation.
Fig. 9. Comparison of cooling performance.
where Tw represents the temperature of the wall near the fluid, Tf represents the temperature of the fluid, and Tc the temperature of the electronic component. According to equations (13) and (14), the heat transfer coefficient h can be expressed as:
1 (Te − Tf ) δ = ∗A− h λ Φ
(15)
In this case, the parameters A, δ and λ were constant. The normal heat transfer coefficient Nh was defined as
Nh =
Φ
(Te − Tf )
(16)
Therefore, Nh can efficiently describe the heat transfer coefficient, and its value will be proportional to h.
We will also assume:
ΔNhi = Nhi−1 − Nhi where i is the number of the electronic component. Although the value of Nh might not reflect the cooling performance difference clearly, ΔNhi could indicate the difference between adjacent electronic components more effectively. An increase in ΔNhi corresponds to a decrease in the cooling performance with respect to the previous component. Fig. 9a and b shows that ΔNh10 was always much larger than ΔNh6 and ΔNh8. In other words, the cooling performance of number 10 decreased sharply. The electronic components named 6 and 8 were both located in the continuous bend channel, whereas number 10 was not. The three components were equipped with the same heat flux, same velocity, same heat transfer area and same geometrical surface, the only possible difference being vapor quality and the
262
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
Δ Nh10
0.3 0.25
1000W
0.2
1600W
0.15 0.1 0.05 0
0
0.005
0.01
Δχ
0.015
0.02
(a)
Fig. 10. Change of ΔNh10 with the change of Δχ.
possible flow behavior. The worst cooling performance of number 10 may either be caused by the large vapor quality or by the different flow behavior. If the vapor quality occupied the main regime of the heat transfer, the decrease of the cooling performance would decrease progressively and ΔNhi would not increase sharply. Since this is not the case, the difference of the cooling performance was most likely caused by the phase distribution. As shown in Fig. 9a and b, the cooling performance of number 10 decreased sharply. The vapor quality under number 8 was smaller than number 10. Fig. 10 shows that the cooling performance difference between 8 and 10 decreases with the increase in Δχ. Therefore, the vapor quality was not the main reason for the sharp decrease of the cooling performance of the electronic component numbered 10. Then, phase distributions have an important influence on heat transfer in a low quality vapor boiling flow. In this case, the flow behavior that was analyzed in the simulation would be confirmed. Actually, a small Δχ corresponded to a weak disturbance of vapor and liquid, and to a large ΔNh10, in agreement with the tendency shown in Fig. 10. It was thus possible to conclude that the numerical analyses were reasonable and supported by the experimental data. 4.3. Optimization of channel and the cooling effect As demonstrated above, the vapor flow tendency might improve the evaporative cooling performance. The effects of low mass flow rate were discussed, due to the fact that such condition can help save energy. Moreover, the pressure loss would increase and the stability of the system would decrease accordingly, if the mass flow was high. Actually, the heat transfer coefficient with a low vapor quality was not as high as the one with a higher quality, according to the results given by Silva Lima et al. [6]. Meanwhile, the low vapor quality burdened the condenser, and it would be possible to save energy if the low mass flow could meet the cooling requirements. Therefore, it was significant to pay more attention to the flow with low mass flow rate. A mass flow rate value of 0.1 t/h was selected for further study. In the coiled channel, the vapor gathered at the inner side of the bend. When the mass flow was low, the vapor velocity was low. If the force could not carry the vapor downstream, the vapor perhaps generated a vortex, leading to a decrease in cooling performance. Since 3-dimensional simulations cost much in terms of computing resources and time [25], further flow analysis was carried out in one bend channel, which was selected in the coiled path. Fig. 11a shows the streamline of the boiling flow. The formation of a vortex near the inner side of the channel was clearly observed.
(b) Fig. 11. Streamline and phase contours of one curve channel.
Comparing with the phase distributions, we can state that most vapors flowed in the area where a vortex formed. In this case, the vapor cannot flow downstream in time, leading to partial retention. Then, the vapor layer and vapor vortex would be in direct contact with the high heat flux surface. Evaporative cooling is known to occur, thanks to the high phase change coefficient. Since the vapor occupied part of the channel and absorbed heat continuously, the cooling performance had to decrease. The forces acting on the vapor include gravity, the pressure force, tension, the Marangoni force, buoyancy and the drag force. Here, the main reason for the generation of a vapor vortex was that the resultant force could not supply enough centrifugal force. If the local velocity of the main streamline was high enough for the drag force to overcome the problem, the vortex would fade away. It was necessary to optimize the model. Hence, two simulations were taken into consideration to confirm whether an increase in local velocity could eliminate the vortex. A contrastive simulation was conducted, which are the original model and an optimal model with a decrease in the effective flow area of the bend (see Table 2). Figs. 12 and 13 show the main streamline and the phase distributions of the original and optimal model, respectively. In the original model, a vortex hid in the main streamline near the inner side of the bend, in agreement with the phase distributions. In the main flow, when the temperature of the fluid increased above the
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
263
Table 2 Working condition of the two models. Model
Mass flux
Saturation temperature
Heat load (W)
Optimal model Original model
0.1 0.1
298 298
100 100
saturation temperature, the vapor emerged and flowed directly downstream, constrained by the above mentioned force. On the one hand, the vapor formed a tendency against the current flow, when the centrifugal force was not enough. On the other hand, it was also driven by the main streamline and pressure force. This led to the vortex generation, and subsequently to partial vapor retention. Complete evaporation in the presence of continuous high flux heating may lead to dangerous effects, and at the same time the cooling performance must worsen. This phenomenon must therefore be dismissed in order to both improve the cooling performance and increase the channel security. In the presence of a vortex, the flow equilibrium can be fully reached, and the velocity balance can be taken by considering all the force acting on the flow. The fluid and vapor flow and the vapor accumulate until the balance is broken. This can also occur when the local velocity is increased with respect to the equilibrium, resulting in the formation of new phase distributions.
(a)
(b) Fig. 13. Streamline and phase contours of optimized model.
(a)
(b) Fig. 12. Streamline and phase contours of original model.
To inquire whether the increase of the local velocity could eliminate the vapor retention, two contrastive models were investigated. The two models were simulated in the same conditions (same mass flow, saturation temperature, pressure, heat load). Fig. 12a shows that the vapor forms a vortex near the inner side of the channel. Five streamlines were selected near the inlet side. In area A, the streamline was basically parallel to the axis of the channel, and the flow was stable, until becoming chaotic when the vapor flowed to area B. When it came to streamline 1 only, the streamline was clearly seen to deviate from the original path. This happened mainly because the forces from the outer side fluid could not supply enough centrifugal force. Then, the streamline bent, developing a vortex in area B. The vapor unfortunately flow in this area (see Fig. 12b), and the high heat flux also transfer through here. The heat transfer performance would worsen and the temperature of the electronic components would increase. Additionally, the streamlines in the outer side area were not parallel to the axis either. Streamline 5 showed a tendency toward the inner side when it flowed through area B. Similarly, the forces, supplied by the fluid in the inner side, were less strong than the opposite force. As a whole, the streamline bent to the inner side of the channel. The open area of the original model near the bend was too large to adapt to the boiling flow here, which directly caused vapor retention near the inner side of the bend, as clearly shown in Fig. 11b. The red color represents the vapor, and the blue shows the liquid. The shape and size were
264
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
similar to the vortex area in the streamline figure, proving that the vortex of the flow field does cause vapor retention. However, Fig. 13a shows a stable flow field whose flow direction is almost parallel to the axis of the channel in area A and area B, which is not in agreement with the original model. Fig. 13b also shows that the vapor flow through the bend was orderly without collection or retention. When the local velocity through area B increased, a new flow balance was established. The forces discussed above supplied the flowing through the bend with the required amount of centrifugal force. Then, the vapor could flow through the bend successfully, avoiding direct contact with the high flux for too long, leading to an improvement in the heat transfer performance. Actually, the shape of the optimal model was based on the streamline of the flow in the original model. When the local velocity increased, the new force was generated and the new force balance was established. The channel size was d (see Fig. 13a) corresponding to the length between streamline 5 and the outer contour of the channel (do in Fig. 12a). This proved that the optimized value d can remove the vortex near the inner side of the tube. We may thus state that the decrease of the open area can improve the stability of the flow and eliminate the vortex that is generated in the vapor area, leading to an improvement in the heat transfer performance. 5. Conclusion A 3-dimensional numerical simulation has been developed to investigate the boiling flow and the phase distribution for the electronic cooling of a moving unit. To account for the movement of the unit, the worst case scenario was chosen. Based on the numerical results, the boiling flow and the phase distributions were investigated. The results showed that the main vapor streamline flowed from the inner side of the one bend to the next one, increasing the turbulence and improving the cooling performance. The phase distributions were influenced by the structure of the channel, gravity, local velocity and other factors. In addition, the numerical analysis and the experimental results indicated that the phase distribution mainly affects the heat transfer in a low vapor quality boiling flow of the coiled channel, not the vapor layer. However, if the forces (gravity, buoyancy, tension, pressure force, or the drag force) cannot supply enough centrifugal force for the vapor, flowing through the bend, a vortex may be generated, decreasing the cooling performance. After an optimized analysis of the shape of the channel, we can conclude that decreasing the effective flow area and increasing the local velocity could eliminate the vortex near the inner side of the bend. Acknowledgements We are grateful to the National Natural Science Foundation of China (51376028) and Program for New Century Excellent Talents in University (NCET-12-0449) for funding this research. Nomenclature r T Sm p F g k SE E Δh v Gf
Vapor radius [mm] Fluid temperature [K] Mass source term Pressure [Pa] Force [N] Gravity acceleration [m·s−2] Curvature of interface Energy source Energy [kJ] Enthalpy difference [kJ] Velocity [m/s] Mass flux per area [kg·s−1·m−2]
M R A h Nh
Relative molecular mass Gas constant Heat transfer area [m2] Heat transfer coefficient [W·m−2·K−1] Normal heat transfer coefficient [W·K−1]
Greek symbols ρ Density, [ kg ⋅ m −3 ] Thermal conductivity [ W ⋅ m −2 ⋅ K −1] λ σ Tension [N·m−1] μ Dynamic viscosity [ Pa ⋅ s] α Volume fraction β Absorbed coefficient γ Phase change coefficient Φ Heat flux [W] δ Thickness d Distance Subscripts sat Saturation condition l Liquid g Gas or vapor total All kinds of force * State near saturation c Cell l-g Liquid translate into vapor g-l Vapor translate into liquid m Mass i Number of electronic component w Wall f Fluid e Electronic component References [1] R. Khodabandeh, Thermal performance of a closed advanced two-phase thermosyphon loop for cooling of radio base stations at different operating conditions, Appl. Therm. Eng. 24 (2004) 2643–2655. [2] A.F. Alia, M.S. El-Genk, Spreaders for immersion nucleate boiling cooling of a computer chip with a central hot spot, Energy Convers. Manag. 53 (2012) 259–267. [3] R. Mongia, K. Masahiro, E. Distefano, et al., Small scale refrigeration system for electronics cooling within a notebook computer. Thermal and Thermomechanical Phenomena in Electronics Systems, Itherm 06, the Tenth Intersociety Conference on San Diego, 751–758, 2006. [4] J.B. Marcinichen, D. Wu, S. Paredes, J.R. Thome, B. Michel, Dynamic flow control and performance comparison of different concepts of two-phase on-chip cooling cycles, Appl. Energy 114 (2014) 179–191. [5] J.B. Marcinichen, J.A. Olivier, V. de Oliveira, J.R. Thome, A review of on-chip micro-evaporation: experimental evaluation of liquid pumping and vapor compression driven cooling systems and control, Appl. Energy 92 (2012) 147–161. [6] R.J. da Silva Lima, J.M. Quibén, J.R. Thome, Flow boiling in horizontal smooth tubes: new heat transfer results for R-134a at three saturation temperatures, Appl. Therm. Eng. 29 (2009) 1289–1298. [7] A. Ma, J. Wei, M. Yuan, J. Fang, Enhanced flow boiling heat transfer of FC-72 on micro-pin-finned surfaces, Heat Mass Transfer 52 (2009) 2925–2931. [8] B. Pulvirenti, A. Matalone, U. Barucca, Boiling heat transfer in narrow channels with offset strip fi ns: application to electronic chipsets cooling, Appl. Therm. Eng. 30 (2010) 2138–2145. [9] B. Kim, B. Sohn, An experimental study of flow boiling in a rectangular channel with offset strip fins, Heat Fluid Flow 27 (2006) 514–521. [10] R. Khodabandeh, Heat transfer in the evaporator of an advanced two-phase thermosyphon loop, Int. J. Refrig. 28 (2005) 190–202. [11] A. Feldman, C. Marvillet, M. Lebouché, Nucleate and convective boiling in plate fin heat exchangers, Heat Mass Transfer 43 (2000) 3433–3442. [12] W. Lee, G. Son, Numerical simulation of boiling enhancement on a microstructured surface, Heat Mass Transfer 38 (2011) 168–173. [13] W. Lee, G. Son, H.Y. Yoon, Direct numerical simulation of flow boiling in a finned microchannel, Heat Mass Transfer 38 (2011) 168–173. [14] R. Zhuan, W. Wang, Simulation of subcooled flow boiling in a micro-channel, Int. J. Refrig. 34 (2011) 781–795. [15] R. Zhuan, W. Wang, Simulation on nucleate boiling in micro-channel, Heat Mass Transfer 53 (2010) 502–512.
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
[16] M. Shoji, Y. Takagi, Bubbling features from a single artificial cavity, Int. J. Heat Mass Transfer 44 (2001) 2763–2776. [17] C.K. Yu, D.C. Lu, T.C. Cheng, Pool boiling heat transfer on artificial micro-cavity surfaces in dielectric fluid FC-72, J. Micromech. Microeng. 16 (2006) 2092–2099. [18] C. Kunkelmann, P. Stephan, Numerical simulation of the transient heat transfer during nucleate boiling of refrigerant HFE-7100, Int. J. Refrig. 33 (2010) 1221–1228. [19] E. Wagner, P. Stephan, O. Koeppen, H. Auracher, High resolution temperature measurements at moving vapor/liquid and vapor/liquid/solid interfaces during bubble growth in nucleate boiling, vol. 260-277, 4th International Berlin Workshop on Transport Phenomena with Moving Boundaries, Berlin, 2007. [20] Y.W. Kuang, W. Wang, R. Zhuan, C.C. Yi, Simulation of boiling flow in evaporator of separate type heat pipe with low heat flux, Ann. Nucl. Energy 75 (2015) 158–167.
265
[21] Z. Yang, X.F. Peng, P. Ye, Numerical and experimental investigation of two phase flow during boiling in a coiled tube, Heat Mass Transfer 51 (2008) 1003–1016. [22] V. Talimi, Y.S. Muzychka, S. Kocabiyik, A review on numerical studies of slug flow hydrodynamics and heat transfer in microtubes and microchannels, Int. J. Multiph. Flow 39 (2012) 88–104. [23] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [24] R.W. Schrage, A Theoretical Study of Interphase Mass Transfer, Columbia University Press, New York, 1953. [25] J.B. Haelssig, A.Y. Tremblay, J. Thibault, S.G. Etemad, Direct numerical simulation of interphase heat and mass transfer in multicomponent vapour–liquid flows, Heat Mass Transfer 53 (2010) 3947–3960.
Contents lists available at ScienceDirect
Applied Thermal Engineering j o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / a p t h e r m e n g
Research Paper
Numerical and experimental investigations of electronic evaporative cooling performance with a coiled channel Xiang Yin a, Feng Cao a,*, Lei Jin a, Bin Hu a, Pengcheng Shu a, Xiaolin Wang b a b
School of Energy and Power Engineering, Xi’an Jiaotong University, 28 Xianning West Road, Xi’an, 710049, China School of Engineering and ICT, University of Tasmania, Private bag 65, Hobart, TAS 7001, Australia
H I G H L I G H T S
• • • • •
Two-phase flow cooling performance of a moving electronic unit was studied. The distributions of the vapor and the liquid phase in a 3D snake pipe were simulated. An experimental setup was established to verify the simulation results. The heat concentration in the curve of the channel affects the cooling effect directly. An optimized method was developed to improve the overall cooling performance.
A R T I C L E
I N F O
Article history: Received 12 July 2015 Accepted 22 October 2015 Available online 3 November 2015 Keywords: Electronic evaporative cooling Numerical simulation Coiled channel VOF Phase distribution Cooling performance
A B S T R A C T
Numerical and experimental investigations on two-phase evaporative cooling performance for a moving electronic unit with high heat flux have been performed in a coiled channel heat exchanger. Based on the VOF model, the phase distribution of a 3D coiled channel was investigated. It was found that the main vapor streamline flowed through the inner side of one bend channel to the next one, which increased the contact frequency of the vapor and the liquid, improving the cooling performance. Moreover, the results were also verified by the experiment. The cooling effect was directly affected by heat concentration in the bend channel. An optimized method was implemented for this phenomenon. Decreasing the effective flow area and increasing the local velocity could establish a new force balance for the vapor and liquid, eliminating the heat concentration and improving the overall cooling performance. © 2015 Elsevier Ltd. All rights reserved.
1. Introduction With the development of human civilization, electronic equipment sets, such as computers, radios, data servers and electric vehicles, intrude all respects of our daily lives. Some particular electronic equipment sets with smaller size have to deal with a large heat flux. Most of them need to work in relatively stable conditions with a small temperature fluctuation range. When the temperature exceeds certain values, the performance of the components goes down sharply, eventually stopping the device from working. Therefore, the cooling performance of electronic components is always an interesting topic for researchers. However, traditional cooling systems can no longer satisfy the requirements for the cooling performance of electronic components because of their limited cooling capability. Nowadays, two-phase evaporative cooling systems are widely applied in the electronic cooling field, e.g. cooling of radio [1], computers [2,3], or chips [4,5]. To improve
* Corresponding author. Tel.: +86 29 82663583; fax: +86 29 82663583. E-mail address: [email protected] (F. Cao). http://dx.doi.org/10.1016/j.applthermaleng.2015.10.127 1359-4311/© 2015 Elsevier Ltd. All rights reserved.
its application, the two-phase boiling performance should be further studied via experiments or numerical simulations. Extensive experimental literature reports on evaporative cooling are available, while only a few deal with evaporative cooling in the macroscopical evaporator by using numerical methods. In a boiling flow, both the convective boiling and the nucleate boiling regimes affect heat transfer. Different mechanisms constitute the main regime in different flows with different shapes [6]. Aixiang Ma et al. [7] concluded that the inner geometry of evaporators has little influences on the two-phase boiling heat transfer when heat loads are high [8]. Moreover, the Reynolds number has a direct effect on heat transfer coefficient of forced convective boiling [9]. In addition, vapor friction has a little influence on the heat transfer coefficient, while the heat flux and pressure have a stronger influence on the heat transfer [10]. That is to say, the effect of inner geometry on the boiling heat transfer is determined by the application and channel structure. A. Feldman et al. [11] studied the mechanism of boiling heat transfer and figured out that the heat transfer coefficient is independent of mass velocity and quality in the convective regime, but is affected by mass velocity and quality in the nucleate regime. In summary, the above mentioned papers
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
mainly showed that the two-phase cooling system could be better employed in the cooling of electronic components. Nevertheless, much research is still needed in view of better applications, especially in larger electronic units. Experimental studies on boiling flow are sometimes limited by practical conditions. It is hard to optimize the structure of the evaporator channel, if the boiling performance cannot meet the requirements. Thus, numerical simulations can be an alternative way to solve the problems. Some efforts in the development of twophase simulations, using the VOF model, have been done to trace the phase interface in a boiling flow. Woorim Lee et al. [12,13] successfully employed the VOF model to analyze the effect of cavity diameter and surface modification on boiling heat transfer. Moreover, numerical simulations have been widely used to investigate the heat transfer performance on the microstructure in a boiling flow [14,15]. The diameter and the number of cavities, as well as the structure of channels, play different roles in different boiling flows [16,17]. The numerical method could be well used to develop the surface configuration for boiling flows. At larger length scales, Christian Kunkelmann and Peter Stephan [18] conducted a numerical investigation on boiling of HFE7100 according to the experiment of Wanger et al [19]. They found that the departure diameter of a bubble is 2.3 mm, and the results are around 20% above the values obtained experimentally. However, theirs was a 2-dimensional simulation, which could not be directly applied into practice. It is worth pointing out that Y.W. Kuang et al. [20] simulated the flow boiling behavior in 3-dimensional pipes, whose diameters were 32 mm and 65 mm, respectively. The model could well simulate the flow pattern evolution, and capture the hydrodynamic as well as thermal mechanisms. But the geometric model is just a straight pipe. In addition, Z. Yang et al. [21] simulated the boiling flow in a coiled tube, and the results were in good agreement with the phenomenon observed in experiments. It is reasonable to claim that in those works the VOF model could be widely used to analyze the phase distribution. However, they did not analyze the correlation between the heat transfer performance and the phase distribution. As a whole, in the field of boiling flow, a number of investigations have been conducted with the VOF model, including the behavior of the bubble, the evolution of the flow pattern, and the performance of heat transfer. However, most of them are limited to a 2-d scale or to an ideal model. Actually, the influence of the bend on flow pattern or heat transfer of two-phase evaporative cooling has been neglected in the past [22]. The numerical investigation of the bend will be highlighted in this paper. In the present paper, an attempt was made to investigate the cooling performance difference of a coiled channel and the vapor behavior, which may affect the cooling performance. A numerical simulation, using the VOF model, and the corresponding experiments were conducted to study the evaporative cooling performance for electronic components in a horizontal coiled heat exchanger with a rectangular channel, which could increase the effective heat transfer area. The heat flux of every component ranges from 60 to 120 W/ cm2. The phase distribution was also studied. Numerical analyses of the cooling performance were verified by experimental data directly and indirectly. Moreover, the channel structure was optimized. 2. Simulation method 2.1. Model of simulation The tracking of phase interface, in the VOF model, was accomplished by solving the continuity equations for the volume fractions of different phases. When the number of the cells is large enough, the errors due to the difference between curves and straight lines could be neglected. The straight line could replace the smooth interface of the two-phase (see Fig. 1). To avoid truncation errors, which would cause a non-conservation of mass velocity between inlet and
257
Fig. 1. Interpolation of interface in VOF.
outlet, the mesh step-size must be selected to match the vapor size. However, in a transient simulation, the time-step is much sensitive to the mesh step-size and must be decreased to confirm stability and convergence, which usually cost much more time. In particular, in 3-dimensional simulations, the simulation period would be too long. Therefore, the mesh step-size should be increased at the expense of simulation accuracy, if the non-conservation is less than 5%. That is to say, some small-sized vapors are neglected.
2.2. Governing equations The fluid in the control volume is governed by three conservation laws: conservation of mass, conservation of momentum and conservation of energy. In the VOF method, the summation of the volume fraction of liquid and vapor equated to 1. The conservation of mass is:
∂α l S + ∇ (v l ⋅ α l ) = m ∂t ρl
(1)
∂α g S + ∇ (v g ⋅ α g ) = m ∂t ρg
(2)
where Sm is the mass source term that reflects the mass transfer between the two phases in a boiling flow, and l, v represent liquid and vapor, respectively. The conservation of momentum is:
∂ ( ρv ) + ∇ ⋅ ( ρvv ) = −∇p + ∇ ⋅ [ μ (∇v + ∇v T )] + ρg + Ftotal ∂t
(3)
The force, Ftotal, is the sum of all forces that affect the generation and movement of vapor, including the surface tension, the force between two phases and so on, excluding gravity. Here, the surface tension is highlighted because of phase changes. In a CSF model [23], the surface tension contributes to surface pressure, generating a force that affects the momentum of the fluid. The force follows that
Fg = σ k
∇ρ ρ ρl − ρ g 1 2 ( ρl + ρ g )
(4)
where k is the curvature of the interface, formed by the liquid and the vapor. It is defined as
k=
Δρ ∇ρ
(5)
258
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
The equation for conservation of energy is:
∂ ( ρE ) + ∇ ⋅ [v ( ρE + p )] = ∇ ⋅ (k eff ∇T ) + S E ∂t
(6)
where SE is the energy source term that satisfies the following equation:
S E = S V ∗ Δh
(7)
where Δh is the enthalpy difference between the liquid and vapor at the state point of saturation temperature in the condition of interface. 2.3. Mass and heat transfer theory
Fig. 2. Geometrical configuration.
In a boiling flow, when the liquid temperature increases above the saturation temperature or the vapor temperature decreases below the saturation temperature, phase change occurs, following mass and heat transfer through the interface of the liquid and vapor to obtain the thermodynamic equilibrium. Based on the kinetic theory, the Hertz Knudsen equation was proposed to calculate mass transfer through the free interface in 1953 [24].
Gf = β
M ( p * − p sat ) 2π RTsat
(8)
where β is the coefficient that reflects the ratio of absorbed vapor to the total amount of vapor entering the liquid surface. The flow was treated as in equilibrium at each time step. The vapor, which is generated from the heat wall, is considered as a constant size globular. Combined with the Clausius Clapeyron equation, the source term per volume (Sm) could be described by the flowing equation:
S m = G f ∗ A c V c = G f ∗ ( πd 2 ) ( 4 πr 3 3) = γ ∗ α g ρg
γ =
( T − Tsat )
⎛ ρl ⎞ 3 M β Δh ⎜ r 2π RTsat ⎝ ρl − ρg ⎟⎠
Tsat
(9)
(10)
3. Experimental facilities
The saturation temperature was treated as a constant number due to the small pressure drop less than 0.1 bar. In a control volume, when the temperature increases above the saturation temperature, some liquid transforms into vapor, to maintain a certain temperature. Conversely, some vapor transforms into liquid to keep the equilibrium. Thus, thermodynamic equilibrium is achieved. The control equation of phase change can be written as follows: If Tc > Tsat, the liquid transforms into vapor:
S m(l − g ) = γ ∗ α l ρl
(T − Tsat ) Tsat
to the fact that such working condition constitutes the worst case in the presence of a rotational movement that changes the position of the cooling equipment. The electronic components were replaced by thermal resistances with the same power that were attached to the double layer heat expanding device in order to decrease the heat flux and improve the cooling performance. Fig. 2 shows the configuration of the electronic components on the cooling heat exchanger with a rectangular coiled channel. The working fluid was R134a, and the secondary boundary was used on the electronic components. To improve the precision and reduce the simulation period, a structured grid was adopted, due to the absence of tortuous surfaces that would require unstructured grids to meet the bodyfitted performance. The mesh step was 0.5 mm in the channel and 0.25 mm near the boundary. The mesh step of the solid flied was 1 mm. In the case of transient simulation, using PISO for pressure interpolation allows lager relaxation factors, decreasing the simulation period. Time step was adapted to the global Courant number to obtain a stable simulation, using values of the order 1E-6 s. The k-ε model was adopted to compute turbulence effects. Additionally, the geometric reconstruction method was used to calculate the phase interface.
To demonstrate the reliability of numerical analysis, experimental investigations were conducted. The experimental equipment is shown in Fig. 3. The mass flow rate of the test heat exchangers was supplied by a sliding-vane working pump. The hand value could increase the adjustable range of the mass flow rate, besides
(11)
Else Tc < Tsat, the vapor transforms into liquid:
S m ( g −l ) = γ ∗ α g ρ g
(T − Tsat ) Tsat
(12)
where γ is calculated according to equation (10). In a certain simulation, γ is constant due to the fact that parameters like ρ, β, and d are definite value. When it comes to R134a, γ was set to 500 or more if the heat flux was too high to maintain thermodynamic equilibrium. 2.4. Geometrical configurations A horizontal working condition was performed to investigate the cooling performance for a moving electronic unit, according
Fig. 3. Experimental facility.
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
adjustment with inverter control. An electric heater, in front of the cooling heat exchanger, could control the vapor quality, which flows into the heat exchanger. The reservoir could improve the stability of the system. Four pressure transducers were placed at the inlet and outlet of the cooling heat exchanger and the pump. Meanwhile, twelve thermocouples were installed near the inlet and outlet of the heat exchanger and every electronic component. 4. Results and discussion 4.1. Phase distribution and cooling performance With the VOF model, the phase distributions could be observed directly, as shown in Fig. 4. The red color represents vapor, and the blue represents liquid. A steady phase distribution was selected to investigate the cooling performance. Although the flow pattern was not steady due to the changing flow parameters (e.g. pressure) in boiling flow, it had little effect on the overall flow tendency when the flow pattern developed well. Therefore, the cooling performance could be successfully investigated from a steady result. The
Fig. 4. Phase distribution of the whole channel.
259
phase distributions present here were selected at a relatively steady flow time. Vapor quality increased along the flow direction, and the amount of liquid decreased. The vapor mainly occupied the top side of the channel, as a consequence of gravity. Originally, the channel was full of liquid. When the electronic components were heated, the heat layer formed. Then, the vapor appeared where the temperature was higher than the saturation temperature. Fig. 4 shows that the vapor tended to gather at the inner side of the bend. Then, in a continuous bend, the main vapor streamline flowed from the inner side of one bend to the inner side of the next one, forming an angle between the main vapor streamline and the channel axis. Meanwhile, the continuous bend made the main vapor streamline change sharply, disturbing the phase interface, which should otherwise be smooth. Fig. 5a and b shows a 3-d translucent main vapor streamline map. The phase interface in the bend is rougher than the one in the straight channel. Thus, heat and mass occur through the whole phase interface form all directions, instead of a single surface. Furthermore, Fig. 5c and e shows the phase interface and the contact action of the liquid and vapor in the continuous bend and straight channel, respectively. Finally, Fig. 5d and f represents the section cut from planes S1 and S2. In the straight channel, the phase interface was smooth under low mass flow condition. Contact between liquid and vapor occurred just below the phase interface. However, in the continuous bend, the rough phase interface and the twisty main vapor streamline improved the contact opportunity of the liquid and the vapor. Mass and heat transfer occurred through the phase interface from all directions, thus increasing the contact frequency of the vapor and the liquid, and improving the cooling performance. Actually, the vapor converging tendency was not clearly observable, because new vapor was continuously generated at the top of the channel, slightly hiding the vapor tendency. Fig. 6a shows the phase contour upon getting rid of the top of the channel, and Fig. 6b represents the sketch map. It is easy to observe that the main vapor streamline does not flow at a constant section. The vapor always gathered at the inner side of the bend. If this vapor could flow downstream in time, the cooling performance would be enhanced; otherwise, a boiling crisis might occur. In our case, the vapor could flow successfully. The cooling performance improved at the continuous bend channel. However, the flow behavior was different at the last bend. The vapor flowed in a smoother manner.
Fig. 5. Difference flow performance in bend and straight channels.
260
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
Table 1 Working conditions of the experiment.
Group 1 Group 2
(a) Phase contour of the continuous
Mass flux/t/h
Total power of components/W
Air temperature/K
0.2~1.6 per o.2 0.2~1.6 per o.2
1000 1600
298 298
tions cut from the direction of SA, SB, and SC in Fig. 4. The main vapor streamline became thin or disappeared, while it was smooth in the channel after the last bend. The flow showed a different behavior, in agreement with the analysis above. Except the last bend, the vapor did not flow at a constant section. The vapor flowed at a constant section of the channel after the last bend, but the phase interface after the last bend was non-turbulent. The contact frequency between the vapor and the liquid was smaller. The cooling performance turned out to be worse than the others. Here, all the flow parameters were the same, except the vapor layer and the phase distribution. The vapor layer might cause a gradual decrease in the cooling performance, but different flow distributions would cause the cooling performance to decrease sharply. Consequently, the cooling performance of the last bend must be much worse than the others. 4.2. Experimental validation
(b) Sketch map of the vapor mainstream Fig. 6. Sketch map of the vapor flow through the continuous curve.
The contact frequency and area at the last bend are expected to be smaller. Therefore, the heat transfer and cooling performances would be worse at this bend than at the others in the same condition (i.e. same velocity, heat flux and heat area, geometrical surface). In addition, it was interesting to find that the vapor flow was not continuous at a constant section. Fig. 7b and c shows the sec-
(a) Phase contour of section SA
(b)
Phase contour of section SB
To demonstrate the analysis above, two groups of experiments were conducted (see Table 1). The electronic components from inlet to outlet were defined as numbers 1, 2, 3…, respectively, and numbers 5, 6, 8, and 10, which had similar heat flux, were selected to analyze the cooling performance. The cooling performance was characterized by the temperature of electronic components: the higher the temperature, the worse cooling performance, and vice versa. 4.2.1. Comparison of experiments and simulations Fig. 8a and b shows the electronic temperature discrepancy between experiment and simulation, whose total heat power is 1000 W and 1600 W, respectively. A similar tendency was found. The electronic component named number 1 might be characterized by a single phase flow, leading to a smaller discrepancy. When the phase changed, the error increased. The error values ranged from 1.23 K to 6.15 K, the largest one taking place under number 3. This may have happened because the subcooled temperature of the simulation is not the same as the experiment, leading to a smaller vapor quality under number 3 than the experiment. The reason why the cooling performance of the simulation was better than the experiment is that the actual thermal contact resistance between the electronic component and the cooling heat exchange was not as ideal as the simulation. The constant globular hypothesis and the simplifications of the VOF model, such as the repetitive computation of the interface advancement, also increased the error. 4.2.2. Indirect validation of the flow pattern The actual flow pattern in the channel was difficult to discover. Here, the flow performance could be deduced from the analysis of the cooling performance. According to the previous analysis, the cooling performance of number 10 was expected to be the worst. The convection heat transfer equation can be written as
Φ = hA ( Tw − Tf
)
(13)
while the heat conduction equation can be written as
(c) Phase contour of section SC Fig. 7. Phase contour of different section.
Φ ( T − Tw ) =λ e δ A
(14)
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
338 336
experiment
0.3
simulation
0.25
ΔNh6
0.2
ΔNh8
0.15
ΔNh10
ΔNh
temperature/K
334 332
0.1
330
0.05
328
0
326 0
2
4
number
6
8
364
(a)
ΔNh6
0.2
Simulation
360
Comparsion of cooling performance under 1000W
0.25
Experiment
362
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 mass flux /t*h-1
10
(a) Temperature comparison under 1000W
ΔNh8 ΔNh10
0.15
358
ΔNh/W*K-1
temperature/K
261
356
0.1
354
0.05
352 350 0
2
4 6 number
8
10
0 0.2
0.4
0.6 0.8 1 1.2 -1 mass flux /t*h
1.4
1.6
(b) Temperature comparison under 1600W (b) Comparison of cooling performance under 1600W
Fig. 8. Temperature comparison of experiment and simulation.
Fig. 9. Comparison of cooling performance.
where Tw represents the temperature of the wall near the fluid, Tf represents the temperature of the fluid, and Tc the temperature of the electronic component. According to equations (13) and (14), the heat transfer coefficient h can be expressed as:
1 (Te − Tf ) δ = ∗A− h λ Φ
(15)
In this case, the parameters A, δ and λ were constant. The normal heat transfer coefficient Nh was defined as
Nh =
Φ
(Te − Tf )
(16)
Therefore, Nh can efficiently describe the heat transfer coefficient, and its value will be proportional to h.
We will also assume:
ΔNhi = Nhi−1 − Nhi where i is the number of the electronic component. Although the value of Nh might not reflect the cooling performance difference clearly, ΔNhi could indicate the difference between adjacent electronic components more effectively. An increase in ΔNhi corresponds to a decrease in the cooling performance with respect to the previous component. Fig. 9a and b shows that ΔNh10 was always much larger than ΔNh6 and ΔNh8. In other words, the cooling performance of number 10 decreased sharply. The electronic components named 6 and 8 were both located in the continuous bend channel, whereas number 10 was not. The three components were equipped with the same heat flux, same velocity, same heat transfer area and same geometrical surface, the only possible difference being vapor quality and the
262
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
Δ Nh10
0.3 0.25
1000W
0.2
1600W
0.15 0.1 0.05 0
0
0.005
0.01
Δχ
0.015
0.02
(a)
Fig. 10. Change of ΔNh10 with the change of Δχ.
possible flow behavior. The worst cooling performance of number 10 may either be caused by the large vapor quality or by the different flow behavior. If the vapor quality occupied the main regime of the heat transfer, the decrease of the cooling performance would decrease progressively and ΔNhi would not increase sharply. Since this is not the case, the difference of the cooling performance was most likely caused by the phase distribution. As shown in Fig. 9a and b, the cooling performance of number 10 decreased sharply. The vapor quality under number 8 was smaller than number 10. Fig. 10 shows that the cooling performance difference between 8 and 10 decreases with the increase in Δχ. Therefore, the vapor quality was not the main reason for the sharp decrease of the cooling performance of the electronic component numbered 10. Then, phase distributions have an important influence on heat transfer in a low quality vapor boiling flow. In this case, the flow behavior that was analyzed in the simulation would be confirmed. Actually, a small Δχ corresponded to a weak disturbance of vapor and liquid, and to a large ΔNh10, in agreement with the tendency shown in Fig. 10. It was thus possible to conclude that the numerical analyses were reasonable and supported by the experimental data. 4.3. Optimization of channel and the cooling effect As demonstrated above, the vapor flow tendency might improve the evaporative cooling performance. The effects of low mass flow rate were discussed, due to the fact that such condition can help save energy. Moreover, the pressure loss would increase and the stability of the system would decrease accordingly, if the mass flow was high. Actually, the heat transfer coefficient with a low vapor quality was not as high as the one with a higher quality, according to the results given by Silva Lima et al. [6]. Meanwhile, the low vapor quality burdened the condenser, and it would be possible to save energy if the low mass flow could meet the cooling requirements. Therefore, it was significant to pay more attention to the flow with low mass flow rate. A mass flow rate value of 0.1 t/h was selected for further study. In the coiled channel, the vapor gathered at the inner side of the bend. When the mass flow was low, the vapor velocity was low. If the force could not carry the vapor downstream, the vapor perhaps generated a vortex, leading to a decrease in cooling performance. Since 3-dimensional simulations cost much in terms of computing resources and time [25], further flow analysis was carried out in one bend channel, which was selected in the coiled path. Fig. 11a shows the streamline of the boiling flow. The formation of a vortex near the inner side of the channel was clearly observed.
(b) Fig. 11. Streamline and phase contours of one curve channel.
Comparing with the phase distributions, we can state that most vapors flowed in the area where a vortex formed. In this case, the vapor cannot flow downstream in time, leading to partial retention. Then, the vapor layer and vapor vortex would be in direct contact with the high heat flux surface. Evaporative cooling is known to occur, thanks to the high phase change coefficient. Since the vapor occupied part of the channel and absorbed heat continuously, the cooling performance had to decrease. The forces acting on the vapor include gravity, the pressure force, tension, the Marangoni force, buoyancy and the drag force. Here, the main reason for the generation of a vapor vortex was that the resultant force could not supply enough centrifugal force. If the local velocity of the main streamline was high enough for the drag force to overcome the problem, the vortex would fade away. It was necessary to optimize the model. Hence, two simulations were taken into consideration to confirm whether an increase in local velocity could eliminate the vortex. A contrastive simulation was conducted, which are the original model and an optimal model with a decrease in the effective flow area of the bend (see Table 2). Figs. 12 and 13 show the main streamline and the phase distributions of the original and optimal model, respectively. In the original model, a vortex hid in the main streamline near the inner side of the bend, in agreement with the phase distributions. In the main flow, when the temperature of the fluid increased above the
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
263
Table 2 Working condition of the two models. Model
Mass flux
Saturation temperature
Heat load (W)
Optimal model Original model
0.1 0.1
298 298
100 100
saturation temperature, the vapor emerged and flowed directly downstream, constrained by the above mentioned force. On the one hand, the vapor formed a tendency against the current flow, when the centrifugal force was not enough. On the other hand, it was also driven by the main streamline and pressure force. This led to the vortex generation, and subsequently to partial vapor retention. Complete evaporation in the presence of continuous high flux heating may lead to dangerous effects, and at the same time the cooling performance must worsen. This phenomenon must therefore be dismissed in order to both improve the cooling performance and increase the channel security. In the presence of a vortex, the flow equilibrium can be fully reached, and the velocity balance can be taken by considering all the force acting on the flow. The fluid and vapor flow and the vapor accumulate until the balance is broken. This can also occur when the local velocity is increased with respect to the equilibrium, resulting in the formation of new phase distributions.
(a)
(b) Fig. 13. Streamline and phase contours of optimized model.
(a)
(b) Fig. 12. Streamline and phase contours of original model.
To inquire whether the increase of the local velocity could eliminate the vapor retention, two contrastive models were investigated. The two models were simulated in the same conditions (same mass flow, saturation temperature, pressure, heat load). Fig. 12a shows that the vapor forms a vortex near the inner side of the channel. Five streamlines were selected near the inlet side. In area A, the streamline was basically parallel to the axis of the channel, and the flow was stable, until becoming chaotic when the vapor flowed to area B. When it came to streamline 1 only, the streamline was clearly seen to deviate from the original path. This happened mainly because the forces from the outer side fluid could not supply enough centrifugal force. Then, the streamline bent, developing a vortex in area B. The vapor unfortunately flow in this area (see Fig. 12b), and the high heat flux also transfer through here. The heat transfer performance would worsen and the temperature of the electronic components would increase. Additionally, the streamlines in the outer side area were not parallel to the axis either. Streamline 5 showed a tendency toward the inner side when it flowed through area B. Similarly, the forces, supplied by the fluid in the inner side, were less strong than the opposite force. As a whole, the streamline bent to the inner side of the channel. The open area of the original model near the bend was too large to adapt to the boiling flow here, which directly caused vapor retention near the inner side of the bend, as clearly shown in Fig. 11b. The red color represents the vapor, and the blue shows the liquid. The shape and size were
264
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
similar to the vortex area in the streamline figure, proving that the vortex of the flow field does cause vapor retention. However, Fig. 13a shows a stable flow field whose flow direction is almost parallel to the axis of the channel in area A and area B, which is not in agreement with the original model. Fig. 13b also shows that the vapor flow through the bend was orderly without collection or retention. When the local velocity through area B increased, a new flow balance was established. The forces discussed above supplied the flowing through the bend with the required amount of centrifugal force. Then, the vapor could flow through the bend successfully, avoiding direct contact with the high flux for too long, leading to an improvement in the heat transfer performance. Actually, the shape of the optimal model was based on the streamline of the flow in the original model. When the local velocity increased, the new force was generated and the new force balance was established. The channel size was d (see Fig. 13a) corresponding to the length between streamline 5 and the outer contour of the channel (do in Fig. 12a). This proved that the optimized value d can remove the vortex near the inner side of the tube. We may thus state that the decrease of the open area can improve the stability of the flow and eliminate the vortex that is generated in the vapor area, leading to an improvement in the heat transfer performance. 5. Conclusion A 3-dimensional numerical simulation has been developed to investigate the boiling flow and the phase distribution for the electronic cooling of a moving unit. To account for the movement of the unit, the worst case scenario was chosen. Based on the numerical results, the boiling flow and the phase distributions were investigated. The results showed that the main vapor streamline flowed from the inner side of the one bend to the next one, increasing the turbulence and improving the cooling performance. The phase distributions were influenced by the structure of the channel, gravity, local velocity and other factors. In addition, the numerical analysis and the experimental results indicated that the phase distribution mainly affects the heat transfer in a low vapor quality boiling flow of the coiled channel, not the vapor layer. However, if the forces (gravity, buoyancy, tension, pressure force, or the drag force) cannot supply enough centrifugal force for the vapor, flowing through the bend, a vortex may be generated, decreasing the cooling performance. After an optimized analysis of the shape of the channel, we can conclude that decreasing the effective flow area and increasing the local velocity could eliminate the vortex near the inner side of the bend. Acknowledgements We are grateful to the National Natural Science Foundation of China (51376028) and Program for New Century Excellent Talents in University (NCET-12-0449) for funding this research. Nomenclature r T Sm p F g k SE E Δh v Gf
Vapor radius [mm] Fluid temperature [K] Mass source term Pressure [Pa] Force [N] Gravity acceleration [m·s−2] Curvature of interface Energy source Energy [kJ] Enthalpy difference [kJ] Velocity [m/s] Mass flux per area [kg·s−1·m−2]
M R A h Nh
Relative molecular mass Gas constant Heat transfer area [m2] Heat transfer coefficient [W·m−2·K−1] Normal heat transfer coefficient [W·K−1]
Greek symbols ρ Density, [ kg ⋅ m −3 ] Thermal conductivity [ W ⋅ m −2 ⋅ K −1] λ σ Tension [N·m−1] μ Dynamic viscosity [ Pa ⋅ s] α Volume fraction β Absorbed coefficient γ Phase change coefficient Φ Heat flux [W] δ Thickness d Distance Subscripts sat Saturation condition l Liquid g Gas or vapor total All kinds of force * State near saturation c Cell l-g Liquid translate into vapor g-l Vapor translate into liquid m Mass i Number of electronic component w Wall f Fluid e Electronic component References [1] R. Khodabandeh, Thermal performance of a closed advanced two-phase thermosyphon loop for cooling of radio base stations at different operating conditions, Appl. Therm. Eng. 24 (2004) 2643–2655. [2] A.F. Alia, M.S. El-Genk, Spreaders for immersion nucleate boiling cooling of a computer chip with a central hot spot, Energy Convers. Manag. 53 (2012) 259–267. [3] R. Mongia, K. Masahiro, E. Distefano, et al., Small scale refrigeration system for electronics cooling within a notebook computer. Thermal and Thermomechanical Phenomena in Electronics Systems, Itherm 06, the Tenth Intersociety Conference on San Diego, 751–758, 2006. [4] J.B. Marcinichen, D. Wu, S. Paredes, J.R. Thome, B. Michel, Dynamic flow control and performance comparison of different concepts of two-phase on-chip cooling cycles, Appl. Energy 114 (2014) 179–191. [5] J.B. Marcinichen, J.A. Olivier, V. de Oliveira, J.R. Thome, A review of on-chip micro-evaporation: experimental evaluation of liquid pumping and vapor compression driven cooling systems and control, Appl. Energy 92 (2012) 147–161. [6] R.J. da Silva Lima, J.M. Quibén, J.R. Thome, Flow boiling in horizontal smooth tubes: new heat transfer results for R-134a at three saturation temperatures, Appl. Therm. Eng. 29 (2009) 1289–1298. [7] A. Ma, J. Wei, M. Yuan, J. Fang, Enhanced flow boiling heat transfer of FC-72 on micro-pin-finned surfaces, Heat Mass Transfer 52 (2009) 2925–2931. [8] B. Pulvirenti, A. Matalone, U. Barucca, Boiling heat transfer in narrow channels with offset strip fi ns: application to electronic chipsets cooling, Appl. Therm. Eng. 30 (2010) 2138–2145. [9] B. Kim, B. Sohn, An experimental study of flow boiling in a rectangular channel with offset strip fins, Heat Fluid Flow 27 (2006) 514–521. [10] R. Khodabandeh, Heat transfer in the evaporator of an advanced two-phase thermosyphon loop, Int. J. Refrig. 28 (2005) 190–202. [11] A. Feldman, C. Marvillet, M. Lebouché, Nucleate and convective boiling in plate fin heat exchangers, Heat Mass Transfer 43 (2000) 3433–3442. [12] W. Lee, G. Son, Numerical simulation of boiling enhancement on a microstructured surface, Heat Mass Transfer 38 (2011) 168–173. [13] W. Lee, G. Son, H.Y. Yoon, Direct numerical simulation of flow boiling in a finned microchannel, Heat Mass Transfer 38 (2011) 168–173. [14] R. Zhuan, W. Wang, Simulation of subcooled flow boiling in a micro-channel, Int. J. Refrig. 34 (2011) 781–795. [15] R. Zhuan, W. Wang, Simulation on nucleate boiling in micro-channel, Heat Mass Transfer 53 (2010) 502–512.
X. Yin et al./Applied Thermal Engineering 94 (2016) 256–265
[16] M. Shoji, Y. Takagi, Bubbling features from a single artificial cavity, Int. J. Heat Mass Transfer 44 (2001) 2763–2776. [17] C.K. Yu, D.C. Lu, T.C. Cheng, Pool boiling heat transfer on artificial micro-cavity surfaces in dielectric fluid FC-72, J. Micromech. Microeng. 16 (2006) 2092–2099. [18] C. Kunkelmann, P. Stephan, Numerical simulation of the transient heat transfer during nucleate boiling of refrigerant HFE-7100, Int. J. Refrig. 33 (2010) 1221–1228. [19] E. Wagner, P. Stephan, O. Koeppen, H. Auracher, High resolution temperature measurements at moving vapor/liquid and vapor/liquid/solid interfaces during bubble growth in nucleate boiling, vol. 260-277, 4th International Berlin Workshop on Transport Phenomena with Moving Boundaries, Berlin, 2007. [20] Y.W. Kuang, W. Wang, R. Zhuan, C.C. Yi, Simulation of boiling flow in evaporator of separate type heat pipe with low heat flux, Ann. Nucl. Energy 75 (2015) 158–167.
265
[21] Z. Yang, X.F. Peng, P. Ye, Numerical and experimental investigation of two phase flow during boiling in a coiled tube, Heat Mass Transfer 51 (2008) 1003–1016. [22] V. Talimi, Y.S. Muzychka, S. Kocabiyik, A review on numerical studies of slug flow hydrodynamics and heat transfer in microtubes and microchannels, Int. J. Multiph. Flow 39 (2012) 88–104. [23] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [24] R.W. Schrage, A Theoretical Study of Interphase Mass Transfer, Columbia University Press, New York, 1953. [25] J.B. Haelssig, A.Y. Tremblay, J. Thibault, S.G. Etemad, Direct numerical simulation of interphase heat and mass transfer in multicomponent vapour–liquid flows, Heat Mass Transfer 53 (2010) 3947–3960.