Nucleate boiling in two types of vertical narrow channels

Abstract

To explore the mechanism of boiling bubble dynamics in narrow channels, we investigate 2-mm wide I- and Z-shaped channels. The influence of wall contact angle on bubble generation and growth is studied using numerical simulation. The relationships between different channel shapes and the pressure drop are also examined, taking into account the effects of gravity, surface tension, and wall adhesion. The wall contact angle imposes considerable influence over the morphology of bubbles. The smaller the wall contact angle, the rounder the bubbles, and the less time the bubbles take to depart from the wall. Otherwise, the bubbles experience more difficulty in departure. Variations in the contact angle also affect the heat transfer coefficient. The greater the wall contact angle, the larger the bubble-covered area. Therefore, wall thermal resistance increases, bubble nucleation is suppressed, and the heat transfer coefficient is lowered. The role of surface tension in boiling heat transfer is considerably more important than that of gravity in narrow channels. The generation of bubbles dramatically disturbs the boundary layer, and the bubble bottom micro-layer can enhance heat transfer. The heat transfer coefficient of Z-shaped channels is larger than that of the I-shaped type, and the pressure drop of the former is clearly higher.

Keywords

nucleate boiling / narrow channel / numerical simulation / bubble dynamics

Introduction

Klausner et al. [1] studied bubble dynamics on walls under flow boiling, measured bubble parameters with a high-speed camera, and established a force model for a bubble. Wambsganss et al. [2] studied the boiling heat transfer of R-113 in mini-channels, and found that nucleate and convective boiling were extremely important to heat transfer. Bowers and Mudawar [3] studied two types of mini-channels with hydraulic diameters (d) of 0.51and 2.54 mm, and found that the heat transfer characteristic of flow boiling in micro/mini-channels differed from the conventional scale. De and Xin [4] first studied the heat transfer characteristics of ethanol in a three-dimensional inner micro-fin heat pipe, and found that the inner micro-fin structure can effectively enhance both boiling and condensing heat transfer. Many of the micro-channel two-phase flow studies also showed that surface tension was an important factor affecting the micro-channel flow pattern [5,6].

Highly compact micro-electronic devices require greater chip cooling capacity, which accelerates the study of the boiling heat transfer mechanism of vapor-liquid two-phase flow in micro-channels. Thus far, however, the heat transfer enhancement mechanism of micro-channel systems has not yet been revealed. In this paper, numerical simulation is used to study the mini-channel boiling phenomenon. The experimental conditions in Ref. [7] are adopted as the prototype using the Fluent software and UDF programming methods to explore the characteristics of boiling bubble dynamics in the two-dimensional mini-channel model. The simulation results are compared with experimental findings.

Numerical simulation process

Mathematical model and solution

Two types of narrow channels, I- and Z-shaped channels, with 2 mm width and 100 mm length, are investigated. The left and right walls are taken as the heating surface. Different wall contact angles are set to explore the influence of the contact angle on bubble dynamics.

The numerical simulation of water boiling heat transfer is intended to explore the course of bubble generation, growth, and departure. Thus, choosing the VOF model is a necessary step in the simulation. The track of phase interface is achieved by solving the equations of a two-phase volume ratio. The Geometric Reconstruction option is chosen, and Body Force is considered. Because Fluent does not contain the model required for boiling calculation, user-defined UDF functions are needed. The UDF program of boiling heat transfer mainly consists of mass transformation from liquid to gas and gas to liquid, latent heat, and variation of surface tension along with temperature. In addition, the fact that the bubble forms only on the heating wall in the UDF must be clearly defined because the actual boiling process is non-homogeneous.

The exchange relationship of vapor-liquid in the UDF is as follows [8]:

When T\ge T_{\mathrm{sat}} (boiling)

R_{\mathrm{L}}=-\lambda {\alpha }_{\mathrm{L}}{\rho }_{\mathrm{L}}{\displaystyle \frac{\left|T-T_{\mathrm{s}\mathrm{a}\mathrm{t}}\right|}{T_{\mathrm{s}\mathrm{a}\mathrm{t}}}},R_{\mathrm{v}}=\lambda {\alpha }_{\mathrm{L}}{\rho }_{\mathrm{L}}{\displaystyle \frac{\left|T-T_{\mathrm{s}\mathrm{a}\mathrm{t}}\right|}{T_{\mathrm{s}\mathrm{a}\mathrm{t}}}}.

When (condensing)

R_{\mathrm{L}}=\lambda {\alpha }_{\mathrm{v}}{\rho }_{\mathrm{v}}{\displaystyle \frac{\left|T-T_{\mathrm{s}\mathrm{a}\mathrm{t}}\right|}{T_{\mathrm{s}\mathrm{a}\mathrm{t}}}},R_{\mathrm{v}}=-\lambda {\alpha }_{\mathrm{v}}{\rho }_{\mathrm{v}}{\displaystyle \frac{\left|T-T_{\mathrm{s}\mathrm{a}\mathrm{t}}\right|}{T_{\mathrm{s}\mathrm{a}\mathrm{t}}}},

where {\alpha }_{\mathrm{v}} and {\alpha }_{\mathrm{L}} are, respectively, volume fractions of the vapor and liquid phases, {\alpha }_{\mathrm{v}}+{\alpha }_{\mathrm{L}}=\mathrm{1}; {\rho }_{\mathrm{v}} and {\rho }_{\mathrm{L}} are, respectively, densities of the vapor and liquid phases, kg/m³; R_{\mathrm{v}} and R_{\mathrm{L}} represent, respectively, source terms of vapor and liquid, kg/(m³·s); T is system temperature, K; T_{\mathrm{s}\mathrm{a}\mathrm{t}} is saturated temperature, K; and \lambda denotes relaxation parameter, s^-1. The heat flux on the interface is q=h_{\mathrm{L}\mathrm{v}}{A}_{\mathrm{L}\mathrm{v}}(T_{\mathrm{v}}-T_{\mathrm{L}})=R_{\mathrm{L}}\times h_{\mathrm{f}\mathrm{g}}, which can be obtained by C-VOF and C-R in Fluent.

Additionally, according to the results of previous experiments [9], the relationship of surface tension and temperature are set as

\sigma =\mathrm{0.09537}-\mathrm{2.24}\times {\mathrm{10}}^{-\mathrm{6}}T-\mathrm{2.560}\times {\mathrm{10}}^{-\mathrm{7}}{T}^{\mathrm{2}}.

Setting up mesh and boundary conditions

The boiling model in a narrow channel is highly sensitive to grids. The thin layer close to the heating wall is the bubble-generating region, and the subtle changes in grid size are likely to lead to non-convergent solutions. In this paper, the thin layer close to the wall undergoes a mesh refinement process, and the minimum mesh size is determined via the critical radius of the bubble, whose radius formula [10] is

r_{\mathrm{c}}={\displaystyle \frac{\mathrm{2}\sigma T_{\mathrm{s}\mathrm{a}\mathrm{t}}}{h_{\mathrm{f}\mathrm{g}}{\rho }_{\mathrm{v}}{(T_{\mathrm{L}}-T_{\mathrm{s}\mathrm{a}\mathrm{t}})}_{\mathrm{c}}}\mathrm{.}}

When the temperature of the heating surface is determined, the bubble radius can be estimated. Based on experience, a minimum mesh size of 1/10 of r_{\mathrm{c}} is appropriate.

Figure 1, the independent validation of the numerical solution, shows that compared with 1.6×10⁶ and 1.8×10⁶ meshes, the relative errors of both T and Dp are within 5%. Thus, the numerical solution of 1.8×10⁶ meshes can be considered grid-independent.

Surface tension acts on the surface, the difference of which causes changes in contact angle with the wall. The surface tension model of Fluent is the continuous surface force model proposed by Brackbill. As the calculation of the triangular and tetrahedral mesh of the surface tension is less accurate than that of quadrilateral and hexahedral grids, quadrilateral and hexahedral meshes should be adopted in the most important regions affected by surface tension.

When the phase interface is tracked using the VOF model, the surface tension must be set. The effect of surface tension is sometimes even larger than that of gravity on the process of boiling heat transfer in mini-channels. Whether the effect of surface tension is important to numerical simulation depends on two non-dimensional Re and Ca(Capillary number), or \mathit{R}\mathit{e} and We (Weber number).

In water boiling, bubbles are produced on the wall. When the wall, water, and steam are in contact, the contact angle is formed between the bubbles and the wall. In this case, when the VOF and surface tension models are both considered, specifying the contact angle of the wall is necessary. To study the effect of the contact angle on boiling, the contact angles are set to 30°, 60°, 120°, and 160° in the calculations.

According to the experimental condition in Ref. [7], the left and right walls serve as the heating surface with constant heat flux. The mass flow rate from the bottom is 83.6 kg/(m²·s), and the upper outlet is free outflow. In the UDF, the system initial temperature is set to 370 K, the boiling temperature at 373.15 K, the inlet pressure at 125 kPa, and the wall toughness at 6.4 μm.

Solution setting

Bubble generation and motion are transient and time-related behaviors. The geometrical-reconstruct method is adopted in the computation. Water is defined as the main phase, and vapor the second phase. The Body Force Formulation option is selected. By balancing the pressure gradient and momentum equation in the volume force, such treatment partially improves the convergence of solutions. Specified Operating Density is chosen and the vapor density is set under Operating Density to eliminate the accumulation of water static pressure and increase round-off accuracy of the momentum balance. The PRESTO method is adopted in pressure interpolation, and PISO is adopted in velocity pressure coupling, without weakening the stability of the solution. Other solution settings can be found in Ref. [11].

Heat and mass transfer on the phase interface

In this paper, the Knudsen number Kn of the established model is far less than 0.001, and the N-S equation is applicable to the simulated object. Heat transfer and mass transfer appear concurrently on the vapor-liquid interface S(r,t) during boiling. Temperature change is continuous but not smooth. The enthalpy, density, and normal seepage velocity are not continuous on the interface.

The mass conservation equation applies in the phase interface.

For the vapor phase

{\displaystyle \frac{\partial }{\partial t}}({\alpha }_{\mathrm{v}}{\rho }_{\mathrm{v}})+\nabla ·({\alpha }_{\mathrm{v}}{\rho }_{\mathrm{v}}{V}_{\mathrm{v}})=\nabla ·({\rho }_{\mathrm{v}}{D}_{\mathrm{v}}\nabla {\alpha }_{\mathrm{v}})+{\mathit{\Gamma }}_{\mathrm{L}\mathrm{v}}.

For the fluid phase

{\displaystyle \frac{\partial }{\partial t}}({\alpha }_{\mathrm{L}}{\rho }_{\mathrm{L}})+\nabla ·({\alpha }_{\mathrm{L}}{\rho }_{\mathrm{L}}{V}_{\mathrm{L}})=\nabla ·({\rho }_{\mathrm{L}}{D}_{\mathrm{L}}\nabla {\alpha }_{\mathrm{L}})+{\mathit{\Gamma }}_{\mathrm{v}\mathrm{L}},

The velocity of evaporation or condensation is given by Equation [12],

{\mathit{\Gamma }}_{\mathrm{L}\mathrm{v}}={\displaystyle \frac{h_{\mathrm{L}\mathrm{v}}{A}_{\mathrm{L}\mathrm{v}}(T_{\mathrm{s}\mathrm{a}\mathrm{t}}-T_{\mathrm{L}})}{h_{\mathrm{f}\mathrm{g}}}}\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{\Gamma }}_{\mathrm{v}\mathrm{L}}={\displaystyle \frac{h_{\mathrm{L}\mathrm{v}}{A}_{\mathrm{L}\mathrm{v}}(T_{\mathrm{L}}-T_{\mathrm{s}\mathrm{a}\mathrm{t}})}{h_{\mathrm{f}\mathrm{g}}}},

where h_{\mathrm{L}\mathrm{v}} is the interface heat transfer coefficient, A_{\mathrm{L}\mathrm{v}} is the interface area of cell volume, T_{\mathrm{s}\mathrm{a}\mathrm{t}} represents the saturation temperature, h_{\mathrm{f}\mathrm{g}} denotes evaporation latent heat, and D is the diffusion coefficient equal to the effective viscosity of phase.

The energy conservation equation is also valid on the interface. In this study, the liquid and vapor phases are assumed incompressible, and the reciprocal of pressure to time can be disregarded.

The energy conservation equation of the vapor phase is

{\displaystyle \frac{\partial }{\partial t}}({\alpha }_{\mathrm{v}}{\rho }_{\mathrm{v}}{c}_{p\mathrm{v}}{T}_{\mathrm{v}})+\nabla ·({\alpha }_{\mathrm{v}}{\rho }_{\mathrm{v}}{c}_{p\mathrm{v}}{V}_{\mathrm{v}}{T}_{\mathrm{v}})=\nabla ·(k_{\mathrm{v}}^{\prime }\nabla T_{\mathrm{v}})+q_{\mathrm{L}\mathrm{v}}+{\mathit{\Gamma }}_{\mathrm{L}\mathrm{v}}{h}_{\mathrm{L}}.

The energy conservation equation of the liquid phase is

{\displaystyle \frac{\partial }{\partial t}}({\alpha }_{\mathrm{L}}{\rho }_{\mathrm{L}}{c}_{p\mathrm{L}}{T}_{\mathrm{L}})+\nabla ·({\alpha }_{\mathrm{L}}{\rho }_{\mathrm{L}}{c}_{p\mathrm{L}}{V}_{\mathrm{L}}{T}_{\mathrm{L}})=\nabla ·(k_{\mathrm{L}}^{\prime }\nabla T_{\mathrm{L}})+q_{\mathrm{v}\mathrm{L}}+{\mathrm{\Gamma }}_{\mathrm{v}\mathrm{L}}{h}_{\mathrm{v}}.

where q_{\mathrm{L}\mathrm{v}} and q_{\mathrm{v}\mathrm{L}} represent the required transformation energy of vapor-liquid two-phase, h_{\mathrm{L}} denotes the enthalpy of the liquid phase, and h_{\mathrm{v}} is the enthalpy of the vapor phase.

Bubble dynamics model

The bubble adhering to the wall is subjected to buoyancy force F_B, pressure F_P, main stream push force F_T, surface tension F_S, resistance F_D, and inertia force F_I [11], as illustrated in Fig. 2(a).

Among these, surface tension F_S and resistance F_D prevent the bubble from departing from the wall, whereas inertia force F_I prevents the bubble from growing. These three forces are regarded as negative forces. On the other hand, buoyancy force F_B, pressure F_P, and main stream push force F_T impel the bubble to depart from the wall. These forces are regarded as positive forces. The equilibrium equation is F_S+F_D+F_I=F_B+F_P+F_T. When the bubble is of a small diameter, the resultant force of buoyancy force F_B and main stream push force F_T is small. However, when the bubble diameter expands, the bubble deforms under the composition of F_B and F_T.

Assume that bubble growth is completed on the heating wall. After the bubble has completely departed, the subjected force is F_D+F_I=F_B + F_T, as shown in Fig. 2(b), and as long as \left|F_{\mathrm{I}}\right|>\mathrm{0}, the bubble moves upward with acceleration.

Simulation results and analysis

Effect of wall contact angle on bubble pattern

Figure 3 demonstrates that left wall contact angle \alpha of the Z-shaped channel is set to 120°, and that of the right wall is 30°. The bubbles generated on both walls present a significant difference: the wetting of the left wall is poor so that the bubbles are flat, whereas the bubbles generated on the right wall are relatively regular circles. The large contact angle promotes the growth of the bubbles close to the wall, and the bubbles are larger than those on the right. In this case, the contact angle is bound to affect the growth rate and departure frequency. Figure 3 clearly shows that the departure diameter on the right wall is smaller and the departure frequency and velocity are higher. Within 15 ms, the bubbles on the right wall have completed the entire process from generation to departure. However, the large bubbles on the left wall have not completed this process. Small bubbles appearing near the left wall are carried from the right wall by incoming streams, rather than the bubbles departing from the left wall.

Figure 4 shows that contact angle \alpha of the left wall of the Z-shaped channel is set to 160°, and that of the right wall is 60°. Compared with Fig. 3, the contact angle on the left wall continues to increase, which causes the generated bubbles to adhere to the wall more firmly; they also experience more difficulty in departing from the wall. In 20 ms, the bubbles show no signs of departure. Compared with Fig. 3, the bubbles generated on the right wall no longer resemble a circle. Because the wall contact angle increases, the wetting properties are reduced, yielding larger contact areas, longer bubble growth time, and less frequent departures.

Bubble nucleation and effect on heat transfer

Bubble nucleation is generated in specific regions. Even at very low flow rates within narrow channels, the liquid flow velocity around the nucleation area is high and may reach 0.3 m/s. Such a high speed stimulates the occurrence of disturbances in the temperature boundary layer near the wall. The disturbance area is centered on the bubble nucleus with a diameter twice that of the bubble diameter. The disturbed liquid rises with horizontal movement, rotation and vibration; thus, a drift region is created at the upper and lower parts of the bubble bottom. The drift region pushes the superheated liquid into the mainstream, and under the effect of incoming flow, the velocity direction of the superheated liquid changes, as depicted in Fig. 5. Severe disturbance is the main factor that compels an increase in heat transfer coefficient. However, the existence of the bubble bottom micro-layer also offers an important contribution to the increase in heat transfer coefficient. Figure 5 shows that the bubble bottom micro-layer flows into the nucleation region along the wall in a direction opposite to the mainstream, and the heat transferred by the micro-layer liquid is greater than the required latent heat of vaporization during bubble generation. The temperature at the bottom of the bubble is lower than the wall temperature; hence, the bubble absorbs the heat from the wall. Nevertheless, the temperature at the top of the bubble is higher than the temperature of the surrounding liquid. Therefore, the bubble releases heat to the liquid. The highest point of the heat transfer coefficient shown in Fig. 5 is the bubble bottom micro-layer.

Relationship between pressure, heat transfer coefficient, and heating time

Figure 6 illustrates that the flow pressure drop of the Z-shaped narrow channel increases significantly by 20% or more on average. The pressure drop is sensitive to channel shape. The initial phase is a single-phase flow, and the difference between the pressure drops of the two channels is small. However, when bubbles begin to generate, the difference becomes more apparent, and the pressure drop of the Z-shaped channel considerably increases, whereas the pressure drop in the I-shaped channel rises relatively slowly.

The change in channel shape is bound to cause variations in the heat transfer coefficient, as shown in Fig. 7. The heat transfer coefficient of the Z-shaped channel is significantly higher than that of the I-shaped channel by 13% on average. Meanwhile, the increase in the heat transfer coefficient enlarges the pressure drop. After the shape changes, the increase in the heat transfer coefficient is less than that of the pressure drop, resulting in poor economy.

We hold that the flow in the narrow channel with a 2-mm width possesses large resistance, and the minor changes in shape would cause a rapid increase in pressure drop. The narrower the channel, the more sensitive the pressure drop is to the shape.

Effect of channel shape on flow heat transfer

Figure 8 presents the flow field distribution in the Z- and I-shaped channels. At the same heating time, the velocity on the surfaces facing the incoming stream changes dramatically in the Z-shaped channel, and appears as an apparent phenomenon of turbulence enhancement, whereas that on other surfaces changes more slowly. The dramatic velocity change on the surfaces facing the incoming stream in the Z-shaped channel considerably enlarges push force F_T, increasing the necessary power departing from the wall, inevitably leading to the increase in departure frequency of the bubbles generated on the wall. The changes in channel shape strengthen the disturbance, yielding an uneven flow field distribution, in contrast to that of the I-shaped channel. The direction of push force F_T is no longer parallel to the wall, which facilitates bubble departure.

Compared with the Z-shaped channel, the flow field distribution of the I-shaped channel is more regular. Although the generation of the bubble on the wall will exert some influence on the flow field distribution, the velocity field remains uniform. For the I-shaped channel, the frequency of bubble generation and departure is lower, and the maximum diameter of the departure bubble is bigger. The decrease in the frequency of bubble generation and departure causes failure to sufficiently strengthen the disturbance of the boundary layer, which is one of the reasons why the heat transfer coefficient of the I-shaped channel is lower than that of the Z-shaped channel.

Comparison with experimental values

Figure 9 shows the comparison between the results in this paper and existing experimental data [5, 6], with a mean error less than 12%. Under different dryness conditions, the results of numerical simulation are slightly higher than the experimental data probably because of the idealized treatment to the flow, assumed as laminar flow. In addition, we define the surface tension between the water and the bubbles as fixed values, without taking the viscous dissipation into account. From Fig. 9, we can observe that the Z-shaped channel has better heat transfer enhancement characteristics.

Conclusions

The following conclusions are drawn:

1) Wall contact angle has a significant effect on bubble shape. The smaller the contact angle, the rounder the bubbles. When the contact angle enlarges, the contact area between the bubble and the wall increases, and the growth rate and bubble departure frequency decrease.

2) The existence of a bubble bottom micro-layer accelerates the convection heat transfer near the wall. The drift region aggravates the disturbance, destroys the thermal boundary layer, and provides a more substantial contribution to the improvement of the boiling heat transfer coefficient.

3) The Z-shaped narrow channel significantly enhances heat transfer. However, comprehensively analyzing the relationship between pressure drop increase and heat transfer enhancement from the perspective of both technology and economy is necessary.

4) The numerical simulation adopts some idealized assumptions such as: the wall contact angle changes linearly and surface tension is a function of temperature, and so on. The difference in contact angles has a significant effect on the process of bubble generation, growth, and departure. Setting contact angles as they actually are is relatively difficult in the simulation process, and more precise experimental data and theoretical analysis are required to pin this process down.

References

[1] Klausner J F, Mei R, Bernhard D M, Zeng L Z. Vapor bubble departure in forced convection boiling. International Journal of Heat and Mass Transfer, 1993, 36(3): 651–662

[2] Wambsganss M W, France D M, Jendrzejczyk J A, Tran T N. Boiling heat transfer in a small diameter tube. ASME Journal of Heat Transfer, 1992, 115(4): 963–972

[3] Bowers M B, Mudawar I. High flux boiling in low flow rate, low pressure drop mini-channel and micro-channel heat sinks. International Journal of Heat and Mass Transfer, 1994, 37(2): 321–332

[4] De J, Xin D M. Heat transfer performance of ethanol three-dimensional inner micro-fin heat pipe. Journal of Chongqing University, 2003, 26(4): 39–41

[5] Mishima K, Hibiki T. Some characteristics of air-water two-phase flow in small diameter vertical tubes. International Journal of Multiphase Flow, 1996, 22(4): 703–712

[6] Chen E F, Li Y Z, Wang S M. Numerical simulation of subcooled boiling water in vertical concentric ammulus under low pressure. Journal of Xi’an Jiaotong University, 2008, 42(7): 855–859

[7] Chen D Q, Pan L M, Yuan D W, Deng J W, Wang X J. Experimental study on bubble growth in vertical narrow channel. Nuclear Power Engineering, 2008, 29(5): 52–55

[8] Talebi S, Abbasi F, Davilu H. A 2D numerical simulation of sub-cooled flow boiling at low-pressure and low-flow rates. Nuclear Engineering and Design, 2009, 239(1): 140–146

[9] Tang C Y. The relationship between surface tension and temperature of water. Journal of Anqing Normal University, 2000, 6(1): 73–74

[10] Lin R T. Boiling Heat Transfer. Beijing: Science Press, 1988

[11] Guo L, Zhang S S, Chen Y Q, Cheng L. Numerical simulation of boiling heat transfer of water in vertical rectangular mini-channels. Applied Mechanics and Materials, 2010, 29-32: 61–67

[12] Hahne E, Grigull U. Boiling Heat Transfer. Beijing: National Defense Industry Press, 1988

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag Berlin Heidelberg