- Email: [email protected]
Simulation on nucleate boiling in micro-channel
International Journal of Heat and Mass Transfer 53 (2010) 502–512
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Simulation on nucleate boiling in micro-channel Rui Zhuan, Wen Wang * Institute of Refrigeration and Cryogenics, School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
a r t i c l e
i n f o
Article history: Received 10 May 2008 Received in revised form 15 January 2009 Accepted 15 August 2009 Available online 12 October 2009 Keywords: Micro-channel Nucleate boiling Marangoni heat transfer Bubble growth Simulation VOF
a b s t r a c t Using the VOF multiphase flow model, numerical simulations are conducted to investigate the nucleate boiling of water in micro-channels. The Marangoni heat transfer through the bubble surface is analyzed, and is compared with the incipient heat flux at the onset of nucleate boiling in micro-channels. The bubble growth in the channel is divided into two stages. At the initial stage, bubble growth is controlled by surface tension, while at the second stage the incipient heat transfer dominated the boiling process. In the results, the full process of bubble generating, growing, departing, combining, and shrinking in the channel is displayed. The simulated results with similar condition are agreed well with some experimental results in references. The method and discussion in the paper are helpful to the investigation of the mechanism of micro-scale two-phase flow and heat transfer. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Two-phase flow and heat transfer in micro-channel has attracted many attentions in recent years. Although there are still some intractable problems, it has been taken into much account because of the perfect performance of heat dissipation and the application of more and more micro-channel heat exchangers. The two-phase flow in a micro-channel is more complicated than the single phase flow, and the nucleate boiling process coupled with high heat flux affects the whole performance significantly. Moreover, in contrast to common scale flow, the surface tension plays a more pronounced function than the gravity, and the stratified two-phase flow was rarely observed. On the other hand, bubbly flow also seldom appears because bubbles grow, coalesce, and reach the channel size very quickly. Thus, some conclusions from the macro-scale can hardly be extrapolated to the micro-scale directly. The intricate characteristics in micro-channels leave us many challenges still [1–3]. In boiling regime, once the Onset of Nucleation Boiling (ONB) occurs in micro-channels, there are some apparent impacts on heat transfer and flow pressure drop. The ONB is also related to the subsequent behaviors of flow boiling, such as the Onset of Significant Void (OSV) and the Departure from Nucleate Boiling (DNB) [4–6]. In some references about ONB, the convective heat transfer was mostly based on Hsu’s model [7] which was developed from pool boiling without considering the convective flow. Some subsequent
* Corresponding author. Tel./fax: +86 21 34206096. E-mail address: [email protected] (W. Wang). 0017-9310/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2009.08.019
models with graphical or numerical procedures derive the ONB criterion [6,8–11], in which, Liu et al. [11] proposed an analytical model from phase equilibrium theory based on the modified Clausius-Clapeyron equation for predicting the ONB of forced convective flow in micro-channels. Bubble growth and departure from the heated wall are significantly affected by the thermal boundary layer near the wall, and its behavior is quite different from that in the common tube since there is high shear stress in the flow direction in a micro-channel. It was believed that the bubble behavior is controlled under heat transfer in its growing stages [12]. It was known from the previous researches [11,13,14] that the incipient heat flux affected the bubbles’ incipient radii and the bubbles’ growth rates. In addition, it is clear that the bubble size is limited by the channel size at least in the cross sections. Moreover, the bubbles’ growth rates properly decrease with the increase in flow Reynolds number, which is the same as in mini-channels [2]. In addition, the buoyancy force is negligible for the bubble departure while the drag force and surface tension are comparatively significant [13]. According to the previous research [15], the size of bubble departure is mainly influenced by drag force and surface tension in micro-channels. About micro-scale two-phase flow, there left some arguments for identifying the macro-to- micro-scale threshold. Several proposals had been put forward for the transition based purely on the hydraulic diameter for non-circular channels [15,16], which were not combined sufficiently with other factors and their function mechanism. In two-phase flow, the bubble growth would be confined by the flowing space in micro-channel. Therefore, the bubble departure size reflects the function of the geometry on
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
503
Nomenclature A area (m2) thermal diffusivity (m2/s) al d bubble diameter (lm) hydraulic diameter (m) De E energy (J/kg) relative error in pressure jump E (DP) volumetric force (N/m3) Fvol volumetric surface tension force (N/m3) F st v g gravity acceleration, 9.81 m/s2 Hc micro-channel height (lm) latent heat (J/kg) hlv h, Dx, Dy, Dz grid size (lm) I third order identity matrix Ja Jakob number k thermal conductivity (W/m K) Ma Marangoni number m mass flux (kg/m3 s) ~ n normal direction to the surface p pressure (Pa) Q heat source term due to phase change (W/m3) wall heat flux (W/cm2) q00 qmarangoni, qm Marangoni heat flux on interface (W/cm2) wall heat flux of onset of nucleate boiling (W/cm2) qONB Re Reynolds number r, R bubble radius (lm) S mass source term due to phase change (kg/m3 s) s tangential direction to the surface t time (s) the previous time in unsteady iteration (s) t0 T temperature (°C)
the flow to some extent. In previous works, the maximum confined single bubble size could also be taken as the micro-scale threshold [1]. Kew and Cornwell [17] recommended a confinement number Co ([r/(g(q‘ qv))]1/2/De) as the criterion to differentiate between macro-scale and micro-scale on two-phase flow and heat transfer with the suggestion neglecting the influence of the flow velocity on the bubble departure diameter. An improved model for departure diameter considering the influence of shear effects has been presented [15]. According to the model, the bubble size increases with the decrease in uout, and there exists the minimum uout for each bubble location. If uout is less than the value, the bubble would grow beyond the confines of the micro-channel, and its shape would no longer be treated as a truncated sphere. As a result, it is reasonable to consider the departure diameter as the macroto-micro-scale threshold. Up to now, the numerical simulation is increasingly and widely employed to investigate two-phase flow. But, in contrast to a macro-scale case, it has been difficult to obtain satisfied simulations on micro-scale flow because of the complexity of the problems and the shortage of feasible models. Still, this work is very necessary since it is attributed to displaying the details in the process of micro-scale boiling. For taking cognizance of nucleate boiling in micro-channels deeply, it is necessary to understand the heat and mass transfer through the bubble surface. Unfortunately, in the previous experimental and theoretical studies, it was not well understood [18]. Thome [1] ever researched the thin film evaporation around elongated bubble and developed a flow boiling model in micro-channels. This model was suitable for elongated bubbles, but not well for the nucleation bubbles in the onset of nucleate boiling in micro-channels.
V wc ww DP
volume (m3) micro-channel width (lm) micro-channel fin thickness (lm) pressure jump (Pa)
Greek symbols volume fraction of liquid phase volume fraction of vapor phase volume fraction of liquid phase (i=l), or vapor phase (i=v) b micro channel aspect ratio ~ t velocity vector (m/s) l dynamic viscosity (kg/ms) q density (kg/m3) r surface tension coefficient j surface curvature g fin efficiency
a‘ av ai
Subscripts b bubble f fluid ‘ liquid min minimum value ONB onset of nucleate boiling sat saturation v vapor vol volume w wall sup superheated degree
At the onset of nucleate boiling, bubbles growth is not confined by channel in general, and the surface tension prevents the bubbles from growing. Thus, the Marangoni flow induced by surface tension would be very important in nucleate boiling [12,19,20]. Arlabosse et al. [21] proved that the Marangoni flow at the nucleation point could intensify the heat transfer on the wall. Moreover, Petrovic et al. [18] evaluated the expression through pool boiling experiment and compared the qm with qONB and divided nucleate boiling process into two stages. In the first stage, it belongs to subcooled boiling while the nucleating bubble on the wall is suppressed by the surface tension. Here, the Marangoni heat flux qm is higher than incipient heat flux qONB, and bubbles growth is controlled by Marangoni heat flow. In the second stage, when the superheat degree of the wall temperature increases constantly and the induced heat flux is large enough for the nucleating bubble to overcome the suppression of surface tension, bubbles are able to depart from the wall and nucleation boiling occurs. That is, in the second stage, the bubble behavior is controlled by incipient heat flux. In simulations of liquid-vapor flow, the volume-of-fluid (VOF) method has been adopted on macro and mini scale boiling flow due to its effect on stability and its efficiency in computation. In order to realize the interface tracking of two-phase flow, the Piecewise Linear Interface Calculation (PLIC) for VOF has been developed [22–28] with consideration of the surface tension force, which can be modeled well with the Continuum Surface Force (CSF) method. Moreover, Brackbill et al. [29,30] recommended addition of a density scaling factor in order to reduce the spurious currents. Recently, some improved methods and algorithm have been carried out by many researchers to increase the accuracy of interface tracking [31–34].
504
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
Bubble behavior in micro-channel is influenced by many factors, such as drag force, surface tension force, lift-off diameter and channel dimension, they need to be considered to investigate the bubble behavior in micro-scale boiling and flowing. The objective of this study is to develop a numerical method to explore the nucleate boiling process in micro-channels. 2. The model for interface tracking on the onset of nucleate boiling The sum of all phase volume fractions should keep a constant,
a‘ þ av ¼ 1
ð1Þ
where al and av are the volume fraction of the liquid phase and vapor phase, respectively. The continuity equations for the volume fractions of all phases are:
@ a‘ S ta‘ Þ ¼ þ r ð~ @t q‘ @ av S tav Þ ¼ þ r ð~ @t qv
ð2Þ ð3Þ
in which S is the source term including the mass transfer rates through bubbles. The mixed density is averaged by the volume fraction, and so are the other mixed properties (conductivity, viscosity, etc.).
q ¼ a‘ q‘ þ av qv
ð4Þ
Momentum equation,
@ tÞ þ r ðq~ t~ tÞ ¼ rp þ r ðq~ @t þ q~ g þ F vol
2 3
lðr~ t þ r~ tT Þ lr ~ tI
3. Heat and mass transfer on the bubble surface in nucleate boiling With the influence of the phase interface, micro bubbles exhibit some unique dynamic behavior [35–37]. In a micro-channel, the interfacial effects on a bubble are usually far more important than the gravitational effects. For example, the buoyancy-to-surface tension ratio is on the order 104 on a water bubble of 100 lm diameter. The bubble was usually supposed to be as a spherical shape with a uniform saturation temperature for simplifying the calculation [38,39]. The bubble is assumed contacting with the wall due to a strong local Marangoni effect, furthermore, the bottom of the bubble is assumed in thermal equilibrium with the wall. In fact, there is a thin liquid film and thermal gradient, which magnitude depends on the temperature profile of the liquid near the wall, sometimes this effect from the gradient on the bubble growth is apparent. Moreover, the vapor and liquid phase are in equilibrium under saturated condition and the heat transfer between the bubble/liquid interfaces is related to the bubble growth rates. In the vicinity of the channel surface, there is an assumed linear temperature gradient profile in the fluid. 3.1. Wall superheat equations in micro-channels In most subcooled boiling, the fluid temperature is lower than the saturation temperature, while the wall temperature is higher than the saturation temperature. The nucleation will occur at some tiny cavities where the local wall temperature is high enough. The nucleating bubble does not depart from the wall until the onset of nucleate boiling (ONB) occurs [18,20]. The wall temperature is dominated by the incipient heat flux [11,40].
pffiffiffiffiffiffiffiffi 2rC b wc þ ww 00 q þ 2 T sat qv h‘v kf 1 þ 2gb Hc sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2r C b wc þ ww 00 q qv h‘v kf 1 þ 2gb Hc
T w ¼ T sat þ ð5Þ
Energy equation,
@ tðqE þ pÞÞ ¼ r ðkrTÞ þ Q ðqEÞ þ r ð~ @t
ð11Þ
ð6Þ 3.2. Heat transfer on the bubble surface in nucleate boiling
here, the heat transfer rates through the interface are embedded in the source term Q. The surface tension force is taken into account with the continuum surface force (CSF) model, and is reformed into an equivalent body force as:
~ F st v ¼ rjn
ð7Þ
The interface curvature and the surface tension are written as below:
~ n ¼ rai
ð8Þ
jðxÞ ¼ r ~ n
ð9Þ
With the phase change at interface, there is a pressure jump at the interface due to the surface tension force, and the two phases will have different pressure. The pressure inside bubble is expressed as following:
n Pv ¼ P‘ þ rj~
ð10Þ 7
8
In the simulation, time steps are chosen as 10 s and 10 s for the calculations, and then the stability criterion is satisfied and the required accuracy is obtained. The grids is hexahedral, and Dx = Dy = Dz = h. The geometric reconstruction interpolation scheme is used in the VOF solution. The Pressure-Implicit with Splitting of Operators (PISO) is considered as the pressure–velocity coupling method, which allows for keeping the solution stability with high values of under-relaxation factors.
Usually, the Knudsen number of the vapor phase is lower than 0.001 in two-phase flow, that is, the flow can be thought of as continuous flow. Supposing that the liquid–gas flow could be considered as being in quasi-equilibrium, the strain is linearly proportional to the stress. Meanwhile, the velocity slip and the temperature skip are neglected on the boundary. 3.2.1. Marangoni heat transfer around bubble In micro-channels, the velocity and temperature gradients near the wall might be extremely large because the liquid is subcooled near inlet. Therefore, thermal capillary force and the drag force resulting from the liquid velocity gradient are significant [41]. Marangoni flow always occurs around the bubble where a temperature gradient exists, meanwhile, the bubble is suppressed by surface tension and attached to the wall surface. Based on equations in Refs. [18,21], Marangoni heat flux can be calculated as following:
qm ¼ ð1 þ 0:00841Ma1=2 Þ
k‘ ðT w T 1 Þ Rb
ð12Þ
where
Ma ¼
ð@ r=@TÞð@T=@sÞRb
l‘ a‘
2
d
ð13Þ
Since the bubble is sticking on the wall in this stage, for simplifying the calculation, the temperature gradient along bubble surface is assumed as the wall temperature gradient,
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
@T s @T l ¼ l @y wall @s Rb
505
ð14Þ
Fig. 1 shows Marangoni heat flux through interface during bubble growth under various superheat degrees near the wall. Here, the incipient bubble diameter is calculated based on Ref. [11], and the superheat degree is defined from Eq. (11). Marangoni heat flux reaches the peak at the beginning of bubble nucleation, and then declines as the bubble grows, as shown in Fig. 1. This can be illustrated that the wall temperature gradient and surface tension gradient reach the maximum at the incipient nucleation point on the wall, which cause the greater Marangoni flow. When superheat layer becomes uniform and thick on the wall, the temperature gradient around a bubble decreases, and so does the Marangoni flux, then the bubble expands gradually. 3.2.2. Analysis the Marangoni heat transfer through bubble surface To analyze the effect of surface tension, the interface Marangoni heat flux is selected to compare with the incipient heat flux under various superheat degrees of fluid near the wall. It is known that Marangoni heat flux decreases with the expansion of bubble in the superheat liquid near the wall. To analyze the incipient heat flux on nucleation point, here, we substitute the incipient bubble radius db with d for simplifying the computation, then,
Ma ¼
ð@ r=@TÞð@T=@sÞRb
l‘ a‘
2
db ¼
ð@ r=@TÞð@T=@yÞwall 2 db l‘ a‘
Fig. 2. Comparison of Marangoni heat flux with qONB at initial nucleation in subcooled boiling. (See above-mentioned references for further information.)
ð15Þ
With the constant inlet velocity, the computed Marangoni flux on incipient nucleation point and incipient wall heat flux are expressed in Fig. 2. In the first stage where the superheat is less than about 8 °C, qmarangoni is greater than qONB. That is, although the wall temperature is higher than the saturation temperature, the bubble growth is suppressed by the surface tension. In the second stage where the superheat is higher than about 8 °C, qONB is higher than qmarangoni, that is, a bubble can overcome the restriction from surface tension, grows rapidly and departs from the wall. In terms of Fig. 2, it is proposed that the transition from the first stage to the second stage occurs when the incipient heat flux reaches a high value, and the value of superheat degree for the transition is around 8 °C in this calculation. Fig. 3 shows the impact of inlet fluid temperature on Marangoni heat flux in subcooled boiling. While the subcooled liquid touches the hot wall at the entrance, the occurrence of Marangoni flow is possible. Increasing the inlet temperature of fluid implies decreas-
Fig. 3. Effect of inlet temperature on Marangoni heat transfer in subcooled boiling.
ing the subcooled degree, which would weaken the Manrangoni heat flux in subcooled boiling. 3.2.3. Model for heat transfer through the phase interface At the first stage, the heat transfer on the interface is dominated by the Marangoni heat transfer, and the heat flux is defined as:
qm ¼ Dm‘—v h‘v
ð16Þ
The second stage is in the nucleate boiling. Assuming bubbles grow in uniformly superheated liquid [38], the interfacial heat transfer, expressed as Eq. (17), is calculated from the bubble growth rate.
q ¼ Dm‘—v h‘v ¼ f ðRÞ h‘v
ð17Þ
3.3. Mass flow through bubble interface in nucleate boiling In the first stage, the interface mass transfer rates can be derived from the heat flux:
Dm‘—v ¼ Fig. 1. Marangoni heat flux through the bubble interface with various superheat degrees near the wall.
qm n Ab 3 qm ¼ n V b h‘v 4 h‘v rb
where n is the number of the nucleation bubbles.
ð18Þ
506
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
At the second stage, bubbles grow rapidly. Their growth rates are expressed as following [38,39].
rffiffiffiffiffiffiffiffiffiffiffiffi 12al t
R ¼ Ja
h = 2 lm in refined region, the number of cells is 1331,817, 2593,674, respectively. The relative deviations of the interface tracking were also evaluated.
ð19Þ
p
4.1. Simulating results for the first stage of subcooled boiling
where al is the thermal diffusivity, and Ja is the Jakob number,
Ja ¼
q‘ cpl ðT w T sat Þ qv h‘v
ð20Þ
0:899
Supposing that the bubble shape was spherical, the mass transfer rates through the bubble surface can be obtained:
Dmb ¼ mb;t0 þ
From Fig. 1, considering that the superheat is 6 °C, the nonlinear curve fitting can be got as following:
qm ¼ 74:173 db
ð24Þ 0:8974
qm;t0 þDt ¼ 83:058 ðdb;t0 þ DdÞ
¼ 83:058 ðdb;t0 þ f ðqm;t0 ÞÞ
dmb dt mb;t0 dt
ð25Þ
ð21Þ
mb ¼ q‘ V b ¼ q‘ f ðRÞ
ð22Þ
dmb dV b df ðRÞ ¼ q‘ ¼ q‘ dt dt dt
ð23Þ
4. The simulation on the nucleate boiling of parallel microchannels Let it be supposed that the flow of water is under the pressure of 1 atm, and the fluid compressibility is neglected. The mean relative roughness in micro-channel surface is less than 1%. Fig 4 shows the shape of the micro-channel of simulation. In order to compare with experimental results in Ref. [11], its dimension is 275 lm (width) 636 lm (height) 25,400 lm (length) as the reference also. The top of the channel is adiabatic, and the heat flux is distributed uniformly at the bottom of the channel. The simulated flow is supposed to be fully developed. With two uniform grids of h = 15 lm, h = 20 lm in the computational domain, the number of cells is 1036,298, 576,584, respectively. Furthermore, the computational domain is arranged with locally refined grids. With the uniform grids of h = 5 lm,
Fig. 4. The micro-channel for simulation.
0:8974
As mentioned above, the incipient heat flux qONB is lower than the Marangoni heat flux qm in the first stage, and the surface tension induces Marangoni flow and predominates the heat transfer on the interface.Fig. 5 shows the flow behavior of the fluid heated in the initial 50 ms, at that time the fluid temperature is lower than the wall temperature. There is a great temperature gradient on the wall, and some hot points appear locally on the wall. Fig. 6 displays the character of bubbles nucleation boiling and flow at the first stage. In general, the bubble nucleates and grows on the wall, and sticks on the wall in the stage. Their behavior is suppressed significantly by surface tension. It can be observed that the phenomenon resulted from simulation is consistent well with the experiment in Ref. [11]. In Fig. 6(a), at the time of t = 850 ms, several nucleating bubbles appear on the wall except at the entrance. Large bubbles nucleate on the bottom surface. The highest value of the local Reynolds number appears in large bubbles, with a certain distance away from the entrance, as shown in Fig. 6(c). In Fig. 6(b), the local velocity vectors also show three modes of bubble flow: (1) the nucleate bubbles stick on the wall; (2) some small bubbles depart and float in the flowing liquid; (3) the large bubbles flow at high speed. In the first mode, small bubbles stick on the wall. The drag force of liquid is negligible, and bubbles are mainly controlled by the surface tension. In the second mode the velocity of bubble is almost the same as the liquid velocity. At this time, with the increase in bubble radius, drag force is large enough to balance the surface tension. Thus bubbles depart and flow with liquid. In the third mode, inertial force increases with bubble radius. Due to the low density and viscosity in vapor, the large inertial force causes the acceleration in bubbles. Therefore, there are velocity differences between large bubbles and the liquid, as shown in Fig. 6(b) and (c), the advent of large bubble is always accompanied by large velocity. Furthermore, this can be illustrated based on the action of pressure gradient in flow field. That is, because of the light density of vapor, the pressure gradient causes high velocity at a large bubble. According to the simulation, the bubble growth is controlled by Mrangoni heat flux resulted from the large temperature gradient near the wall. Since the heat flux from the wall is not high enough in the calculation, the bubble grows slowly, and no bubble departs
Fig. 5. The temperature (K) field near the entrance for q00 = 16.16 W/cm2, t = 50 ms. h = 15 lm, Rmin/h = 2.
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
507
Fig. 6. Flowing parameters near the entrance at the first stage (qincipient = 16.16 W/cm2, qincipient is lower than qmarangoni.). h = 15 lm, Rmin/h = =2. (a) Volume fraction of bubbles over bottom surface (t = 850 ms), (b) velocity vectors (t = 850 ms), (c) cell Reynolds number (t = 850 ms).
from the wall. The liquid nearby the wall is heated gradually to the nucleation temperature. In addition, there is no spherical bubble with the diameter larger than 200 lm on the wall, since the superheated zone near the wall is very thin. 4.2. The simulation for nucleate boiling in the second stage While the bubble occurs on the hot wall, the heat transfer resulted from the Marangoni effect is helpful to the growth of micro bubbles [43,44]. Within a certain range, the size of the bubble is dominated by the heat flux. In the second stage, the incipient heat flux is higher than the Marangoni heat flux, and then the bubbles grow rapidly. In Fig. 2, when the superheat degree near the wall is over about 8 °C, the Marangoni heat flux at the nucleation point is lower than the incipient heat flux.
Fig. 7. Comparison of Marangoni heat flux in subcooled boiling as well as incipient heat flux.
In Fig. 7, with high incipient heat flux above the transitional point, the fluid temperature increases along the flowing direction and reaches the saturation temperature near the exit, at the same time, the temperature difference near the wall decreases gradually, and the Marangoni heat flux falls below the incipient heat flux. Meantime some bubbles overcome the suppression of surface tension and lift-off from the wall. In Figs. 8 and 9, the heat flux on the wall is q00 = 70 W/cm2, and the liquid is subcooled at entrance. Near the exit, the fluid reaches saturation temperature. It is a typical boiling process, in Fig. 8, the temperature field in the channel is not uniform due to the influence of successive bubbles. Fig. 9 shows the volume fraction variation in the channel, the bubbles are easy to be distinguished from the liquid. In the channel, the whole process of the behavior of bubble is displayed, nucleating on the wall near the entrance, bubble growing rapidly, then bubble lifting off from the wall, and some bubbles are combined into big bubbles. The process of the bubbles lift-off and coalescence at the second stage is shown in Fig. 10, which comes from the vertical symmetric section of the channel. During the initial process (t = 0–200 ms), several bubbles expand, lift-off and combine with each other because of the high value of superheat near the wall. During the process from t = 600 ms to t = 1200 ms, more bubbles coalesce into large bubbles, although the shape of the large bubble is distorted at the beginning. At the time of t = 1700 ms, the large bubble is less distorted due to the curvature correction of the interface tracking. In this figure, it lasts about 1100 ms that the distorted large bubble turn to approximate round shape. In the second stage, under the high wall heat flux, the liquid temperature increases and its temperature gradient near the wall decreases. Thus the effect of Marangoni flow abates. From the simulating results, it can be seen that bubble grows more quickly than in the first stage. While the bubble nucleates on the wall, the wall temperature at this location rises in a similar manner as in the macro-scale boiling.
Fig. 8. The temperature field near the entrance (K) (q00 = 70 W/cm2). h = 15 lm, Rmin/h = 2.
508
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
Fig. 9. Bubbles volume fraction on the bottom surface for q00 = 70 W/cm2. h = 15 lm, Rmin/h = 2.
It is believed that the heat dissipation path from the channel wall to liquid is partially blocked by the vapor bubble. Therefore, the vapor temperature rises with the increase of bubble size. 4.3. Bubbles growth rates Figs. 10 and 11 show bubbles growth and departure in microchannel. At the initial growth process, inertial force is negligible and the shape of micro bubbles is round, as shown in Fig. 11. However, inertial force increases with a bubble radius. As the inertial
force is large enough to balance the surface tension at interface, large bubbles deform and break, as shown in Fig. 10. Fig. 12 shows the simulation about bubbles growth with heating density of 70 W/cm2 and 16.16 W/cm2 during the initial process of nucleation. It can be observed that the incipient heat flux influences the bubble growth significantly. At low incipient heat flux (qincipient = 16.16 W/cm2), the surface tension controls the bubble growth. At this time Marangoni heat flow induced by temperature gradient dominates the bubble interface heat transfer. At high incipient heat flux (70 W/cm2), the surface tension hardly impact the bubbles. Meanwhile the Marangoni flow is weakened due to the little temperature gradient, and the bubble growth rates controlled by high incipient heat flux are kept high values. There were some published experimental data in similar condition in Ref. [11], and the results are also plotted in the diagram, which are consistent very well with the simulated data. Fig. 12 shows the simulation results are influenced by grid size. When bubble radius is higher than 25 lm, the difference between the results increases with calculation time. 4.4. Analysis of errors In the computation, the error in the pressure jump is investigated to estimate the accuracy of the interface tracking. The relative pressure jump error is evaluated as following:
EðDPÞ ¼
jDP DPexact j DPexact
ð26Þ
where DP is the numerical jump in pressure through bubble interface. The exact jump in pressure across a bubble is given by:
DPexact ¼ rj
ð27Þ
With two locally refined grids of h = 5 lm (Rmin/h = 2), h = 2 lm (Rmin/h = 5) in computational region, the number of total grids is 1331,817, 2593,674 respectively, and errors are plotted in Fig. 13(a). With the uniform grids of h = 15 lm (Rmin/h = 2), Fig. 10. Bubbles departure and coalescence at the second stage on the vertical symmetric profiles for q00 = 70 W/cm2. h = 20 lm, Rmin/h = 2.
Fig. 11. Axisymmetric and round bubbles (R < 40 lm) for q00 = 25 W/cm2. h = 5 lm, Rmin/h = 2.
Fig. 12. Comparing the growth of bubble radius resulted from simulation with the experimental data in Ref. [10].
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
509
Fig. 14. Volume fraction of vapor at the outlet of the channel (h = 15 lm, Rmin/h = 2).
16.16 W/cm2), bubble grows slowly, and the volume fraction of vapor at exit is low as well. 4.6. Bubble departure diameter and macro-to-micro transition threshold
Fig. 13. Relative error in calculation of pressure during bubbles growth for q00 = 70 W/cm2. (a) h = 5 lm (Rmin/h = 2), h = 2 lm (Rmin/h = 5); (b) h = 15 lm (Rmin/h = 2), and h = 20 lm (Rmin/h = 2).
h = 20 lm (Rmin/h = 2) in the computational domain, the number of total grids is 1036,298, 576,584 respectively, then errors are plotted in Fig. 13(b). Fig. 13(a) shows the errors in pressure jump, which increase with the grid size and bubble radius. Micro bubble requires much fine grids in order to obtain a solution with reasonable accuracy due to the significant refinement dependency. In addition, when a bubble grows, errors in curvature increase with time. Therefore, errors in pressure jump at interface depend primarily on grid size and time of bubble growth. Moreover, spurious velocities result from errors in curvature estimation. In Fig. 13(b), for a large bubble, errors increase quickly with bubble radius. With the bubble radius increase, drag force of liquid influences the shape of a bubble significantly, and the bubble is deformed easily, as shown in Fig. 10. Therefore, for the large deformed bubbles the refined grids may not increase the accuracy efficiently inside the wobbling region. 4.5. Vapor volume fraction at outlet Fig. 14 shows the variation of volume fraction of vapor at the outlet. The fluctuation comes from the bubble movement in the channel. There is a high value of the mean volume fraction of vapor at the outlet, it is apparent that high heat flux induces the high bubble growth rates. At low incipient heat flux (for example,
In simulation, the flow was supposed to be 1 atm at the entrance of the channel, and the velocity is 0.65 m/s. The departure diameter is 120 lm based on the criteria from Ref. [15], which is less than the channel dimension. In Figs. 6 and 10, there is no apparently jamming with one single bubble on the channel cross section. Although some small bubbles have collected into several large bubbles, their shapes are the truncated sphere still. The bubbles grow gradually due to heating in the flow. Decreasing uout is helpful to their expanding rapidly [15]. There should be the minimum value of uout for each bubble, below the value the bubble could grow beyond the confines of the micro-channel, and its shape would no longer be treated as a truncated sphere. In this simulation, the minimum uout is 0.35 m/ s, which is lower than the mean flowing velocity (0.65 m/s) at the entrance. The bubble behavior takes an important role in flow and heat transfer in micro-channel. It is proposed that the transition criterion of macro-to-micro flow might be related to the bubble departure diameter in two-phase area [1]. That is, while the departure bubble could expand to fill the channel, the further growth would be confined by the channel, the flow should be micro two-phase flow. In this paper, the dbub is less than the size of micro-channels, the shape of the single bubble is a truncated sphere, not the confined elongated bubble. 5. The simulation on the nucleate boiling in a micro trapezoid channel 5.1. Boundary conditions and source term With etching on silicon, the trapezoid channel is very popular, some work on micro-channels are on this kind of channel. For verifying the simulation with above model and method, it is shown that some calculated results are compare with the published data in Ref. [13]. The boundary conditions and fluid properties are quoted from the data in the reference. The micro-channel is 28 mm in length and the cross section is shown in Fig. 15. The
510
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
Fig. 15. The cross section of a single trapezoid micro-channel.
top of the channel is adiabatic, and the uniform heat flux is distributed at the bottom of the channel. The wall temperature can be obtained from Eq. (11). The incipient heat flux q00ONB in experiment [13] is higher than the calculated q00m , so the nucleate boiling should be controlled by the incipient heat flux q00ONB and the heat and mass rates are derived from Eqs. (21) and (23). The computational domain is arranged with locally refined grids, the grid size h is 2 lm (Rmin/h = 4) in refined domain, the number of total grids is 1092,905. 5.2. Results and discussion 5.2.1. The simulation for bubble nucleation and growth Fig. 16 displays the sequential simulated results reflecting the bubbles nucleation and their evolution at the bottom surface at the distance of 2 mm from the outlet. Here, the bubble evolution includes the process of bubble nucleation, growth and departure. Most of the bubbles are not confined from the channel dimension in this stage. Since the dimensions of the channel are different in width and height, the bubble departure diameter is properly greater than the height of the channel, that is, the bubble will be restricted in the vertical dimension in the growing process. Fig. 17 displays the bubble growth and the confined bubble on the vertical section.
The bubbles nucleate from the vicinity of the inlet to the outlet. At isotropic stage bubbles grow with a spherical shape. When the radii of bubbles reach the value of 20–25 lm, they begin to be distorted and become anisotropic due to the limitation of channel depth and the dragging of the liquid. The non-uniform growth rates among bubbles are not surprising as the local flow field may be changing due to the perturbation of the growth of other bubbles at other sites. Fig. 18(a) displays the simulated bubbles radii on position 4, 5, 6 in Fig. 16, which are compared with the experimental data in Ref. [13]. The results for bubbles on position 2, 3 while the heating density q00 is equal to 255 kW/m2 are compared with the experiments in Fig. 18(b). In the experiment, five bubble evolutions were observed on different points for different heat flux. Each process was named as cycle in the figure. The critical radius for the unconfined bubbles is approximately 22–25 lm in terms of simulation, they are closed to the values 25–30 lm from experiments [13]. Furthermore, the three simulated growth rates of bubbles (4)–(6) are the same as each other, and agree well with the referred tested data. The differences among the three positions along the flow are mainly due to the tiny differences on the enthalpy on the points. Fig. 19 shows relative errors in pressure jump. At the initial time of nucleation, bubble shape is controlled by surface tension, and errors increase slightly. When bubble boundary is close to the top wall of channel, due to limit of the channel wall and drag force of liquid, the bubble shape becomes distorted, and then errors increase quickly. It is concluded that when bubbles are disturbed by additional force, relative errors in pressure jump can not be used to estimate the accuracy of the interface tracking. Fig. 20 shows the variation of the calculated average bubbles radii with q00 = 255 kW/m2 and q00 = 189 kW/m2, which are compared with the experimental data in literature. According to the simulation, the heat flux influences the bubbles’ incipient radii and growth rates. However, it is not conspicuous in the experimental results. The experimental data in both heat fluxes correspond to the simulated profile with 255 kW/m2 of heat flux in the growing period. Meanwhile, in the initial nucleating period, the observed bubble radii are larger than the simulated values.
Fig. 16. Bubbles nucleation and growth at the bottom surface downstream for one cycle (t = 0–140 ms), G = 341 kg/m2, q00 = 189 kW/m2, h = 2 lm, Rmin/h = 4.
Fig. 17. Bubbles nucleation and growth at the vertical wall upstream (t = 0–300 ms), G = 341 kg/m2 s, q00 = 189 kW/m2, h = 2 lm, Rmin/h = 4.
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
511
Fig. 20. Time evolution of average bubble radius (h = 2 lm, Rmin/h = 4).
Fig. 18. Time evolution of bubble radius for G = 341 kg/m2 s: (a) q00 = 189 kW/m2, (b) q00 = 255 kW/m2; h = 2 lm, Rmin/h = 4. Fig. 21. Pressure variation at the horizontal section near the outlet (h = 2 lm, Rmin/h = 4).
the maximum pressure that can be sustained in the bubble is dictated significantly by the saturation pressure corresponding to the wall temperature. As shown in Fig. 21, the pressure spike occurs at the nucleation sites which cause the pressure variation along the flow direction. The pressure jump at the nucleation site leads to a disruption in the flow. When bubbles is sticking on the wall, the surface tension force along the contact line may balance the drag force of the incoming liquid, and can disturb flow at interface. 6. Summary
Fig. 19. Relative error in calculation of pressure during bubbles growth for G = 341 kg/m2 s, q00 = 189 kW/m2, h = 2 lm, Rmin/h = 4.
The nucleate boiling in micro-channels was investigated with numerical simulation, and some simulation results were compared with previous experiments in literature. The comparison indicates that the simulation is feasible and reliable. The conclusions are as following:
5.2.2. Flow stability during nucleate boiling in a micro-channel According to Eq. (10), the pressure jump occurs at a bubble interface due to the surface tension. As a bubble grows on a cavity,
(1) Through the analysis with Marangoni heat transfer on bubble surface, it is proposed that Marangoni heat transfer is important at the onset of nucleate boiling in micro-channels. Moreover, the nucleate boiling process is divided into two
512
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
stages. In the first stage, bubbles are suppressed by surface tension, and Marangoni convection dominates the interface heat transfer. In the second stage, the wall superheat temperature is higher than the transitional point, and qONB is large enough to overcome the surface tension. At this time, bubbles grow rapidly and lift off from the wall. (2) The feature of nucleate boiling in common channel is free bubble flow, while in micro-channels it is a confined bubble flow. It is proposed that the macro-to-micro transition criterion in two-phase flow and heat transfer should be related to the bubble departure diameter, furthermore, the width and height might influence the transitional point too. The departure diameter could be considered as the transitional criterion of macro-to-micro dimension. (3) It is proved that the incipient heat flux would affect the bubble incipient radius and the bubble growth rates. In addition, the higher the heat flux, the higher the volume fraction of vapor at the outlet of a micro-channel. (4) The pressure spike occurs at the nucleation sites due to the influence of surface tension on interface. The pressure perturbation at the nucleation sites might be the causes of disruption in the flow.
Acknowledgement The authors appreciate the support from National Natural Science Fund in China (50576054). References [1] John R. Thome, Boiling in microchannels: a review of experiment and theory, Int. J. Heat Fluid Flow 25 (2004) 128–139. [2] S.G. Kandlikar, Fundamental issues related to flow boiling in minichannels and microchannels, Exp. Thermal Fluid Sci. 26 (2002) 389–407. [3] A.E. Bergles, V.J.H. Lienhard, G.E. Kendall, P. Griffith, Boiling and evaporation in small diameter channels, Heat Transfer Eng. 24 (2003) 18–40. [4] S.M. Ghiaasiaan, R.C. Chedester, Boiling incipience in microchannels, Inter. J. Heat Mass Transfer 45 (2002) 4599–4606. [5] R.C. Chedester, S.M. Ghiaasiaan, A proposed mechanism for hydrodynamicallycontrolled onset of significant void in microtubes, Int. J. Heat Fluid Flow 23 (2002) 769–775. [6] J.E. Kennedy, G.M. Roach Jr., M.E. Dowling, S.I. Abdel-Khalik, S.M. Ghiaasiaan, S.M. Jeter, Z.H. Qureshi, The onset of flow instability in uniformly heated horizontal microchannels, ASME J. Heat Transfer 122 (2000) 118–125. [7] Y.Y. Hsu, On the size range of active nucleation cavities on a heating surface, J. Heat Transfer 84 (1962) 207–216. [8] S.G. Kandlikar, V.R. Mizo, M.D. Cartwright, Ikenze, Bubble nucleation and growth characteristics in subcooled flow boiling of water, HTD-vol. 342, in: ASME Proceedings of the 32nd National Heat Transfer Conference, vol. 4, 1997, pp. 11–18. [9] W. Qu, I. Mudawar, Prediction and measurement of incipient boiling heat flux in microchannel heat sinks, Inter. J. Heat Mass Transfer 45 (2002) 3933–3945. [10] J. Li, P. Cheng, Bubble cavitation in a microchannel, Int. J. Heat Mass Transfer 47 (2004) 2689–2698. [11] Dong Liu, Poh-Seng Lee, Suresh V. Garimella, Prediction of the onset of nucleate boiling in microchannel flow, Int. J. Heat Mass Transfer 48 (2005) 5134–5149. [12] N.O. Young, J.S. Goldstein, M.J. Block, The motion of bubbles in a vertical temperature gradient, J. Fluid Mech. 6 (1959) 350–356. [13] P.C. Lee, F.G. Tseng, Chin Pan, Bubble dynamics in microchannels. Part I: single microchannel, Int. J. Heat Mass Transfer 47 (2004) 5575–5589. [14] H.Y. Li, F.G. Tseng, Chin Pan, Bubble dynamics in microchannels. Part II: two parallel microchannels, Int. J. Heat Mass Transfer 47 (2004) 5591–5601. [15] Weilin Qu, Issam Mudawar, Prediction and measurement of incipient boiling heat flux in micro-channel heat sinks, Int. J. Heat Mass Transfer 45 (2002) 3933–3945.
[16] S.S. Mehendal, A.M. Jacobi, R.K. Shah, Fluid flow and heat transfer at micro- and meso-scales with application to heat exchanger design, Appl. Mech. Rev. 53 (7) (2000) 175–193. [17] P. Kew, K. Cornwell, Correlations for prediction of boiling heat transfer in small diameter channels, Appl. Thermal Eng. 17 (8–10) (1997) 705– 715. [18] Sanja Petrovic, Tony Robinson, Ross L. Judd, Marangoni heat transfer in subcooled nucleate pool boiling, Int. J. Heat Mass Transfer 47 (2004) 5115– 5128. [19] J. Li, G.P. Peterson, Microscale heterogeneous boiling on smooth surfaces— from bubble nucleation to bubble dynamics, Int. J. Heat Mass Transfer 48 (2005) 4316–4332. [20] J.L. McGrew, F.L. Bamford, T.R. Rehm, Marangoni flow: an additional mechanism in boiling heat transfer, Science 153 (1966) 1106–1107. [21] P. Arlabosse, L. Tadrist, H. Tadrist, J. Pantaloni, Experimental analysis of the heat transfer induced by thermocapillary convection around a bubble, J. Heat Transfer 122 (2000) 66–73. [22] W.F. Noh, P.R. Woodward, SLIC (Simple Line Interface Method), in: Lecture Notes in Physics, vol. 59, Springer, 1976. [23] C.W. Hirt, B.D. Nichols, Volume of Fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [24] D.L. Youngs, Time-dependent multi-material flow with large fluid distribution, in: K.W. Morton, M.L. Norman (Eds.), Numerical Methods for Fluid Dynamics, 1986, pp. 187–221. [25] N. Ashgriz, J.Y. Poo, FLAIR: flux line-segment model for advection and interface reconstruction, J. Comput. Phys. 93 (1991) 449–468. [26] S.O. Kim, H.C. No, Second order model for free surface convection and interface reconstruction, Int. J. Numer. Methods Fluids 26 (1) (1998) 79–100. [27] J.E. Polliod, An analysis of piecewise linear interface reconstruction algorithm for volume of fluid methods, M.Sc. Thesis, Department of Mathematics, University of California, Davis, 1992. [28] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (1998) 112–152. [29] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [30] J.U. Brackbill, D.B. Kothe, Dynamic modeling of the surface tension, in: Proceedings of the 3rd Microgravity Fluid Physics Conference, Cleveland, OH, 1996, pp. 693–698. [31] M. Seifollahi, E. Shirani, N. Ashgrizb, An improved method for calculation of interface pressure force in PLIC-VOF methods, Eur. J. Mech. B/Fluids 27 (2008) 1–23. [32] E. Shirani, N. Ashgriz, J. Mostaghimi, Interface pressure calculation based on conservation of momentum for front capturing methods, J. Comput. Phys. 203 (2005) 154–175. [33] Daniel Lörstad, Laszlo Fuchs, High-order surface tension VOF-model for 3D bubble flows with high density ratio, J. Comput. Phys. 200 (2004) 153– 176. [34] Marianne M. Francois, Sharen J. Cummins, Edward D. Dendy, Douglas B. Kothe, James M. Sicilian, Matthew W. Williams, A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (2006) 141–173. [35] H. Wang, X.F. Peng, W.K. Lin, C. Pan, B.X. Wang, Bubble-top jet flow on micro wires, Int. J. Heat Mass Transfer 47 (2004) 2891–2900. [36] H. Wang, X.F. Peng, B.X. Wang, D.J. Lee, Bubble-sweeping mechanisms, Sci. China, Ser. E 46 (2003) 225–233. [37] Hao Wang, Xiaofeng Peng, Suresh V. Garimella, David M. Christopher, Microbubble return phenomena during subcooled boiling on small wires, Int. J. Heat Mass Transfer 50 (2007) 163–172. [38] J.R. Thome, V. Dupont, A.M. Jacobi, Heat transfer model for evaporation in microchannels. Part I: presentation of the model, Int. J. Heat Mass Transfer 47 (2004) 3375–3385. [39] Nengli Zhang, David F. Chat, Models for enhanced boiling heat transfer by unusual Marangoni effects under microgravity conditions, Int. Comm. Heat Mass Transfer 26 (1999) 1081–1090. [40] Satish G. Kandlikar, Nucleation characteristics and stability considerations during flow boiling in microchannels, Exp. Therm. Fluid Sci. 30 (2006) 411– 417. [41] C.L. Vandervort, A.E. Bergles, M.K. Jensen, Heat transfer mechanisms in very high heat flux subcooled boiling, ASME Fundam. Subcooled Flow Boiling HTD 217 (1992) 1–9. [42] E.J. Davis, G.H. Anderson, The incipience of nucleate boiling in forced convection flow, AIChE J. 12 (1966) 774–780. [43] Jr. Hung Tsai, Liwei Lin, Transient thermal bubble formation on polysilicon micro-resisters, J. Heat Transfer 124 (2002) 375–382. [44] Liwei Lin, Microscale thermal bubble formation: thermophysical phenomena and applications, Microscale Thermophys. Eng. 2 (1998) 71–85.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Simulation on nucleate boiling in micro-channel Rui Zhuan, Wen Wang * Institute of Refrigeration and Cryogenics, School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China
a r t i c l e
i n f o
Article history: Received 10 May 2008 Received in revised form 15 January 2009 Accepted 15 August 2009 Available online 12 October 2009 Keywords: Micro-channel Nucleate boiling Marangoni heat transfer Bubble growth Simulation VOF
a b s t r a c t Using the VOF multiphase flow model, numerical simulations are conducted to investigate the nucleate boiling of water in micro-channels. The Marangoni heat transfer through the bubble surface is analyzed, and is compared with the incipient heat flux at the onset of nucleate boiling in micro-channels. The bubble growth in the channel is divided into two stages. At the initial stage, bubble growth is controlled by surface tension, while at the second stage the incipient heat transfer dominated the boiling process. In the results, the full process of bubble generating, growing, departing, combining, and shrinking in the channel is displayed. The simulated results with similar condition are agreed well with some experimental results in references. The method and discussion in the paper are helpful to the investigation of the mechanism of micro-scale two-phase flow and heat transfer. Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction Two-phase flow and heat transfer in micro-channel has attracted many attentions in recent years. Although there are still some intractable problems, it has been taken into much account because of the perfect performance of heat dissipation and the application of more and more micro-channel heat exchangers. The two-phase flow in a micro-channel is more complicated than the single phase flow, and the nucleate boiling process coupled with high heat flux affects the whole performance significantly. Moreover, in contrast to common scale flow, the surface tension plays a more pronounced function than the gravity, and the stratified two-phase flow was rarely observed. On the other hand, bubbly flow also seldom appears because bubbles grow, coalesce, and reach the channel size very quickly. Thus, some conclusions from the macro-scale can hardly be extrapolated to the micro-scale directly. The intricate characteristics in micro-channels leave us many challenges still [1–3]. In boiling regime, once the Onset of Nucleation Boiling (ONB) occurs in micro-channels, there are some apparent impacts on heat transfer and flow pressure drop. The ONB is also related to the subsequent behaviors of flow boiling, such as the Onset of Significant Void (OSV) and the Departure from Nucleate Boiling (DNB) [4–6]. In some references about ONB, the convective heat transfer was mostly based on Hsu’s model [7] which was developed from pool boiling without considering the convective flow. Some subsequent
* Corresponding author. Tel./fax: +86 21 34206096. E-mail address: [email protected] (W. Wang). 0017-9310/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2009.08.019
models with graphical or numerical procedures derive the ONB criterion [6,8–11], in which, Liu et al. [11] proposed an analytical model from phase equilibrium theory based on the modified Clausius-Clapeyron equation for predicting the ONB of forced convective flow in micro-channels. Bubble growth and departure from the heated wall are significantly affected by the thermal boundary layer near the wall, and its behavior is quite different from that in the common tube since there is high shear stress in the flow direction in a micro-channel. It was believed that the bubble behavior is controlled under heat transfer in its growing stages [12]. It was known from the previous researches [11,13,14] that the incipient heat flux affected the bubbles’ incipient radii and the bubbles’ growth rates. In addition, it is clear that the bubble size is limited by the channel size at least in the cross sections. Moreover, the bubbles’ growth rates properly decrease with the increase in flow Reynolds number, which is the same as in mini-channels [2]. In addition, the buoyancy force is negligible for the bubble departure while the drag force and surface tension are comparatively significant [13]. According to the previous research [15], the size of bubble departure is mainly influenced by drag force and surface tension in micro-channels. About micro-scale two-phase flow, there left some arguments for identifying the macro-to- micro-scale threshold. Several proposals had been put forward for the transition based purely on the hydraulic diameter for non-circular channels [15,16], which were not combined sufficiently with other factors and their function mechanism. In two-phase flow, the bubble growth would be confined by the flowing space in micro-channel. Therefore, the bubble departure size reflects the function of the geometry on
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
503
Nomenclature A area (m2) thermal diffusivity (m2/s) al d bubble diameter (lm) hydraulic diameter (m) De E energy (J/kg) relative error in pressure jump E (DP) volumetric force (N/m3) Fvol volumetric surface tension force (N/m3) F st v g gravity acceleration, 9.81 m/s2 Hc micro-channel height (lm) latent heat (J/kg) hlv h, Dx, Dy, Dz grid size (lm) I third order identity matrix Ja Jakob number k thermal conductivity (W/m K) Ma Marangoni number m mass flux (kg/m3 s) ~ n normal direction to the surface p pressure (Pa) Q heat source term due to phase change (W/m3) wall heat flux (W/cm2) q00 qmarangoni, qm Marangoni heat flux on interface (W/cm2) wall heat flux of onset of nucleate boiling (W/cm2) qONB Re Reynolds number r, R bubble radius (lm) S mass source term due to phase change (kg/m3 s) s tangential direction to the surface t time (s) the previous time in unsteady iteration (s) t0 T temperature (°C)
the flow to some extent. In previous works, the maximum confined single bubble size could also be taken as the micro-scale threshold [1]. Kew and Cornwell [17] recommended a confinement number Co ([r/(g(q‘ qv))]1/2/De) as the criterion to differentiate between macro-scale and micro-scale on two-phase flow and heat transfer with the suggestion neglecting the influence of the flow velocity on the bubble departure diameter. An improved model for departure diameter considering the influence of shear effects has been presented [15]. According to the model, the bubble size increases with the decrease in uout, and there exists the minimum uout for each bubble location. If uout is less than the value, the bubble would grow beyond the confines of the micro-channel, and its shape would no longer be treated as a truncated sphere. As a result, it is reasonable to consider the departure diameter as the macroto-micro-scale threshold. Up to now, the numerical simulation is increasingly and widely employed to investigate two-phase flow. But, in contrast to a macro-scale case, it has been difficult to obtain satisfied simulations on micro-scale flow because of the complexity of the problems and the shortage of feasible models. Still, this work is very necessary since it is attributed to displaying the details in the process of micro-scale boiling. For taking cognizance of nucleate boiling in micro-channels deeply, it is necessary to understand the heat and mass transfer through the bubble surface. Unfortunately, in the previous experimental and theoretical studies, it was not well understood [18]. Thome [1] ever researched the thin film evaporation around elongated bubble and developed a flow boiling model in micro-channels. This model was suitable for elongated bubbles, but not well for the nucleation bubbles in the onset of nucleate boiling in micro-channels.
V wc ww DP
volume (m3) micro-channel width (lm) micro-channel fin thickness (lm) pressure jump (Pa)
Greek symbols volume fraction of liquid phase volume fraction of vapor phase volume fraction of liquid phase (i=l), or vapor phase (i=v) b micro channel aspect ratio ~ t velocity vector (m/s) l dynamic viscosity (kg/ms) q density (kg/m3) r surface tension coefficient j surface curvature g fin efficiency
a‘ av ai
Subscripts b bubble f fluid ‘ liquid min minimum value ONB onset of nucleate boiling sat saturation v vapor vol volume w wall sup superheated degree
At the onset of nucleate boiling, bubbles growth is not confined by channel in general, and the surface tension prevents the bubbles from growing. Thus, the Marangoni flow induced by surface tension would be very important in nucleate boiling [12,19,20]. Arlabosse et al. [21] proved that the Marangoni flow at the nucleation point could intensify the heat transfer on the wall. Moreover, Petrovic et al. [18] evaluated the expression through pool boiling experiment and compared the qm with qONB and divided nucleate boiling process into two stages. In the first stage, it belongs to subcooled boiling while the nucleating bubble on the wall is suppressed by the surface tension. Here, the Marangoni heat flux qm is higher than incipient heat flux qONB, and bubbles growth is controlled by Marangoni heat flow. In the second stage, when the superheat degree of the wall temperature increases constantly and the induced heat flux is large enough for the nucleating bubble to overcome the suppression of surface tension, bubbles are able to depart from the wall and nucleation boiling occurs. That is, in the second stage, the bubble behavior is controlled by incipient heat flux. In simulations of liquid-vapor flow, the volume-of-fluid (VOF) method has been adopted on macro and mini scale boiling flow due to its effect on stability and its efficiency in computation. In order to realize the interface tracking of two-phase flow, the Piecewise Linear Interface Calculation (PLIC) for VOF has been developed [22–28] with consideration of the surface tension force, which can be modeled well with the Continuum Surface Force (CSF) method. Moreover, Brackbill et al. [29,30] recommended addition of a density scaling factor in order to reduce the spurious currents. Recently, some improved methods and algorithm have been carried out by many researchers to increase the accuracy of interface tracking [31–34].
504
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
Bubble behavior in micro-channel is influenced by many factors, such as drag force, surface tension force, lift-off diameter and channel dimension, they need to be considered to investigate the bubble behavior in micro-scale boiling and flowing. The objective of this study is to develop a numerical method to explore the nucleate boiling process in micro-channels. 2. The model for interface tracking on the onset of nucleate boiling The sum of all phase volume fractions should keep a constant,
a‘ þ av ¼ 1
ð1Þ
where al and av are the volume fraction of the liquid phase and vapor phase, respectively. The continuity equations for the volume fractions of all phases are:
@ a‘ S ta‘ Þ ¼ þ r ð~ @t q‘ @ av S tav Þ ¼ þ r ð~ @t qv
ð2Þ ð3Þ
in which S is the source term including the mass transfer rates through bubbles. The mixed density is averaged by the volume fraction, and so are the other mixed properties (conductivity, viscosity, etc.).
q ¼ a‘ q‘ þ av qv
ð4Þ
Momentum equation,
@ tÞ þ r ðq~ t~ tÞ ¼ rp þ r ðq~ @t þ q~ g þ F vol
2 3
lðr~ t þ r~ tT Þ lr ~ tI
3. Heat and mass transfer on the bubble surface in nucleate boiling With the influence of the phase interface, micro bubbles exhibit some unique dynamic behavior [35–37]. In a micro-channel, the interfacial effects on a bubble are usually far more important than the gravitational effects. For example, the buoyancy-to-surface tension ratio is on the order 104 on a water bubble of 100 lm diameter. The bubble was usually supposed to be as a spherical shape with a uniform saturation temperature for simplifying the calculation [38,39]. The bubble is assumed contacting with the wall due to a strong local Marangoni effect, furthermore, the bottom of the bubble is assumed in thermal equilibrium with the wall. In fact, there is a thin liquid film and thermal gradient, which magnitude depends on the temperature profile of the liquid near the wall, sometimes this effect from the gradient on the bubble growth is apparent. Moreover, the vapor and liquid phase are in equilibrium under saturated condition and the heat transfer between the bubble/liquid interfaces is related to the bubble growth rates. In the vicinity of the channel surface, there is an assumed linear temperature gradient profile in the fluid. 3.1. Wall superheat equations in micro-channels In most subcooled boiling, the fluid temperature is lower than the saturation temperature, while the wall temperature is higher than the saturation temperature. The nucleation will occur at some tiny cavities where the local wall temperature is high enough. The nucleating bubble does not depart from the wall until the onset of nucleate boiling (ONB) occurs [18,20]. The wall temperature is dominated by the incipient heat flux [11,40].
pffiffiffiffiffiffiffiffi 2rC b wc þ ww 00 q þ 2 T sat qv h‘v kf 1 þ 2gb Hc sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2r C b wc þ ww 00 q qv h‘v kf 1 þ 2gb Hc
T w ¼ T sat þ ð5Þ
Energy equation,
@ tðqE þ pÞÞ ¼ r ðkrTÞ þ Q ðqEÞ þ r ð~ @t
ð11Þ
ð6Þ 3.2. Heat transfer on the bubble surface in nucleate boiling
here, the heat transfer rates through the interface are embedded in the source term Q. The surface tension force is taken into account with the continuum surface force (CSF) model, and is reformed into an equivalent body force as:
~ F st v ¼ rjn
ð7Þ
The interface curvature and the surface tension are written as below:
~ n ¼ rai
ð8Þ
jðxÞ ¼ r ~ n
ð9Þ
With the phase change at interface, there is a pressure jump at the interface due to the surface tension force, and the two phases will have different pressure. The pressure inside bubble is expressed as following:
n Pv ¼ P‘ þ rj~
ð10Þ 7
8
In the simulation, time steps are chosen as 10 s and 10 s for the calculations, and then the stability criterion is satisfied and the required accuracy is obtained. The grids is hexahedral, and Dx = Dy = Dz = h. The geometric reconstruction interpolation scheme is used in the VOF solution. The Pressure-Implicit with Splitting of Operators (PISO) is considered as the pressure–velocity coupling method, which allows for keeping the solution stability with high values of under-relaxation factors.
Usually, the Knudsen number of the vapor phase is lower than 0.001 in two-phase flow, that is, the flow can be thought of as continuous flow. Supposing that the liquid–gas flow could be considered as being in quasi-equilibrium, the strain is linearly proportional to the stress. Meanwhile, the velocity slip and the temperature skip are neglected on the boundary. 3.2.1. Marangoni heat transfer around bubble In micro-channels, the velocity and temperature gradients near the wall might be extremely large because the liquid is subcooled near inlet. Therefore, thermal capillary force and the drag force resulting from the liquid velocity gradient are significant [41]. Marangoni flow always occurs around the bubble where a temperature gradient exists, meanwhile, the bubble is suppressed by surface tension and attached to the wall surface. Based on equations in Refs. [18,21], Marangoni heat flux can be calculated as following:
qm ¼ ð1 þ 0:00841Ma1=2 Þ
k‘ ðT w T 1 Þ Rb
ð12Þ
where
Ma ¼
ð@ r=@TÞð@T=@sÞRb
l‘ a‘
2
d
ð13Þ
Since the bubble is sticking on the wall in this stage, for simplifying the calculation, the temperature gradient along bubble surface is assumed as the wall temperature gradient,
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
@T s @T l ¼ l @y wall @s Rb
505
ð14Þ
Fig. 1 shows Marangoni heat flux through interface during bubble growth under various superheat degrees near the wall. Here, the incipient bubble diameter is calculated based on Ref. [11], and the superheat degree is defined from Eq. (11). Marangoni heat flux reaches the peak at the beginning of bubble nucleation, and then declines as the bubble grows, as shown in Fig. 1. This can be illustrated that the wall temperature gradient and surface tension gradient reach the maximum at the incipient nucleation point on the wall, which cause the greater Marangoni flow. When superheat layer becomes uniform and thick on the wall, the temperature gradient around a bubble decreases, and so does the Marangoni flux, then the bubble expands gradually. 3.2.2. Analysis the Marangoni heat transfer through bubble surface To analyze the effect of surface tension, the interface Marangoni heat flux is selected to compare with the incipient heat flux under various superheat degrees of fluid near the wall. It is known that Marangoni heat flux decreases with the expansion of bubble in the superheat liquid near the wall. To analyze the incipient heat flux on nucleation point, here, we substitute the incipient bubble radius db with d for simplifying the computation, then,
Ma ¼
ð@ r=@TÞð@T=@sÞRb
l‘ a‘
2
db ¼
ð@ r=@TÞð@T=@yÞwall 2 db l‘ a‘
Fig. 2. Comparison of Marangoni heat flux with qONB at initial nucleation in subcooled boiling. (See above-mentioned references for further information.)
ð15Þ
With the constant inlet velocity, the computed Marangoni flux on incipient nucleation point and incipient wall heat flux are expressed in Fig. 2. In the first stage where the superheat is less than about 8 °C, qmarangoni is greater than qONB. That is, although the wall temperature is higher than the saturation temperature, the bubble growth is suppressed by the surface tension. In the second stage where the superheat is higher than about 8 °C, qONB is higher than qmarangoni, that is, a bubble can overcome the restriction from surface tension, grows rapidly and departs from the wall. In terms of Fig. 2, it is proposed that the transition from the first stage to the second stage occurs when the incipient heat flux reaches a high value, and the value of superheat degree for the transition is around 8 °C in this calculation. Fig. 3 shows the impact of inlet fluid temperature on Marangoni heat flux in subcooled boiling. While the subcooled liquid touches the hot wall at the entrance, the occurrence of Marangoni flow is possible. Increasing the inlet temperature of fluid implies decreas-
Fig. 3. Effect of inlet temperature on Marangoni heat transfer in subcooled boiling.
ing the subcooled degree, which would weaken the Manrangoni heat flux in subcooled boiling. 3.2.3. Model for heat transfer through the phase interface At the first stage, the heat transfer on the interface is dominated by the Marangoni heat transfer, and the heat flux is defined as:
qm ¼ Dm‘—v h‘v
ð16Þ
The second stage is in the nucleate boiling. Assuming bubbles grow in uniformly superheated liquid [38], the interfacial heat transfer, expressed as Eq. (17), is calculated from the bubble growth rate.
q ¼ Dm‘—v h‘v ¼ f ðRÞ h‘v
ð17Þ
3.3. Mass flow through bubble interface in nucleate boiling In the first stage, the interface mass transfer rates can be derived from the heat flux:
Dm‘—v ¼ Fig. 1. Marangoni heat flux through the bubble interface with various superheat degrees near the wall.
qm n Ab 3 qm ¼ n V b h‘v 4 h‘v rb
where n is the number of the nucleation bubbles.
ð18Þ
506
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
At the second stage, bubbles grow rapidly. Their growth rates are expressed as following [38,39].
rffiffiffiffiffiffiffiffiffiffiffiffi 12al t
R ¼ Ja
h = 2 lm in refined region, the number of cells is 1331,817, 2593,674, respectively. The relative deviations of the interface tracking were also evaluated.
ð19Þ
p
4.1. Simulating results for the first stage of subcooled boiling
where al is the thermal diffusivity, and Ja is the Jakob number,
Ja ¼
q‘ cpl ðT w T sat Þ qv h‘v
ð20Þ
0:899
Supposing that the bubble shape was spherical, the mass transfer rates through the bubble surface can be obtained:
Dmb ¼ mb;t0 þ
From Fig. 1, considering that the superheat is 6 °C, the nonlinear curve fitting can be got as following:
qm ¼ 74:173 db
ð24Þ 0:8974
qm;t0 þDt ¼ 83:058 ðdb;t0 þ DdÞ
¼ 83:058 ðdb;t0 þ f ðqm;t0 ÞÞ
dmb dt mb;t0 dt
ð25Þ
ð21Þ
mb ¼ q‘ V b ¼ q‘ f ðRÞ
ð22Þ
dmb dV b df ðRÞ ¼ q‘ ¼ q‘ dt dt dt
ð23Þ
4. The simulation on the nucleate boiling of parallel microchannels Let it be supposed that the flow of water is under the pressure of 1 atm, and the fluid compressibility is neglected. The mean relative roughness in micro-channel surface is less than 1%. Fig 4 shows the shape of the micro-channel of simulation. In order to compare with experimental results in Ref. [11], its dimension is 275 lm (width) 636 lm (height) 25,400 lm (length) as the reference also. The top of the channel is adiabatic, and the heat flux is distributed uniformly at the bottom of the channel. The simulated flow is supposed to be fully developed. With two uniform grids of h = 15 lm, h = 20 lm in the computational domain, the number of cells is 1036,298, 576,584, respectively. Furthermore, the computational domain is arranged with locally refined grids. With the uniform grids of h = 5 lm,
Fig. 4. The micro-channel for simulation.
0:8974
As mentioned above, the incipient heat flux qONB is lower than the Marangoni heat flux qm in the first stage, and the surface tension induces Marangoni flow and predominates the heat transfer on the interface.Fig. 5 shows the flow behavior of the fluid heated in the initial 50 ms, at that time the fluid temperature is lower than the wall temperature. There is a great temperature gradient on the wall, and some hot points appear locally on the wall. Fig. 6 displays the character of bubbles nucleation boiling and flow at the first stage. In general, the bubble nucleates and grows on the wall, and sticks on the wall in the stage. Their behavior is suppressed significantly by surface tension. It can be observed that the phenomenon resulted from simulation is consistent well with the experiment in Ref. [11]. In Fig. 6(a), at the time of t = 850 ms, several nucleating bubbles appear on the wall except at the entrance. Large bubbles nucleate on the bottom surface. The highest value of the local Reynolds number appears in large bubbles, with a certain distance away from the entrance, as shown in Fig. 6(c). In Fig. 6(b), the local velocity vectors also show three modes of bubble flow: (1) the nucleate bubbles stick on the wall; (2) some small bubbles depart and float in the flowing liquid; (3) the large bubbles flow at high speed. In the first mode, small bubbles stick on the wall. The drag force of liquid is negligible, and bubbles are mainly controlled by the surface tension. In the second mode the velocity of bubble is almost the same as the liquid velocity. At this time, with the increase in bubble radius, drag force is large enough to balance the surface tension. Thus bubbles depart and flow with liquid. In the third mode, inertial force increases with bubble radius. Due to the low density and viscosity in vapor, the large inertial force causes the acceleration in bubbles. Therefore, there are velocity differences between large bubbles and the liquid, as shown in Fig. 6(b) and (c), the advent of large bubble is always accompanied by large velocity. Furthermore, this can be illustrated based on the action of pressure gradient in flow field. That is, because of the light density of vapor, the pressure gradient causes high velocity at a large bubble. According to the simulation, the bubble growth is controlled by Mrangoni heat flux resulted from the large temperature gradient near the wall. Since the heat flux from the wall is not high enough in the calculation, the bubble grows slowly, and no bubble departs
Fig. 5. The temperature (K) field near the entrance for q00 = 16.16 W/cm2, t = 50 ms. h = 15 lm, Rmin/h = 2.
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
507
Fig. 6. Flowing parameters near the entrance at the first stage (qincipient = 16.16 W/cm2, qincipient is lower than qmarangoni.). h = 15 lm, Rmin/h = =2. (a) Volume fraction of bubbles over bottom surface (t = 850 ms), (b) velocity vectors (t = 850 ms), (c) cell Reynolds number (t = 850 ms).
from the wall. The liquid nearby the wall is heated gradually to the nucleation temperature. In addition, there is no spherical bubble with the diameter larger than 200 lm on the wall, since the superheated zone near the wall is very thin. 4.2. The simulation for nucleate boiling in the second stage While the bubble occurs on the hot wall, the heat transfer resulted from the Marangoni effect is helpful to the growth of micro bubbles [43,44]. Within a certain range, the size of the bubble is dominated by the heat flux. In the second stage, the incipient heat flux is higher than the Marangoni heat flux, and then the bubbles grow rapidly. In Fig. 2, when the superheat degree near the wall is over about 8 °C, the Marangoni heat flux at the nucleation point is lower than the incipient heat flux.
Fig. 7. Comparison of Marangoni heat flux in subcooled boiling as well as incipient heat flux.
In Fig. 7, with high incipient heat flux above the transitional point, the fluid temperature increases along the flowing direction and reaches the saturation temperature near the exit, at the same time, the temperature difference near the wall decreases gradually, and the Marangoni heat flux falls below the incipient heat flux. Meantime some bubbles overcome the suppression of surface tension and lift-off from the wall. In Figs. 8 and 9, the heat flux on the wall is q00 = 70 W/cm2, and the liquid is subcooled at entrance. Near the exit, the fluid reaches saturation temperature. It is a typical boiling process, in Fig. 8, the temperature field in the channel is not uniform due to the influence of successive bubbles. Fig. 9 shows the volume fraction variation in the channel, the bubbles are easy to be distinguished from the liquid. In the channel, the whole process of the behavior of bubble is displayed, nucleating on the wall near the entrance, bubble growing rapidly, then bubble lifting off from the wall, and some bubbles are combined into big bubbles. The process of the bubbles lift-off and coalescence at the second stage is shown in Fig. 10, which comes from the vertical symmetric section of the channel. During the initial process (t = 0–200 ms), several bubbles expand, lift-off and combine with each other because of the high value of superheat near the wall. During the process from t = 600 ms to t = 1200 ms, more bubbles coalesce into large bubbles, although the shape of the large bubble is distorted at the beginning. At the time of t = 1700 ms, the large bubble is less distorted due to the curvature correction of the interface tracking. In this figure, it lasts about 1100 ms that the distorted large bubble turn to approximate round shape. In the second stage, under the high wall heat flux, the liquid temperature increases and its temperature gradient near the wall decreases. Thus the effect of Marangoni flow abates. From the simulating results, it can be seen that bubble grows more quickly than in the first stage. While the bubble nucleates on the wall, the wall temperature at this location rises in a similar manner as in the macro-scale boiling.
Fig. 8. The temperature field near the entrance (K) (q00 = 70 W/cm2). h = 15 lm, Rmin/h = 2.
508
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
Fig. 9. Bubbles volume fraction on the bottom surface for q00 = 70 W/cm2. h = 15 lm, Rmin/h = 2.
It is believed that the heat dissipation path from the channel wall to liquid is partially blocked by the vapor bubble. Therefore, the vapor temperature rises with the increase of bubble size. 4.3. Bubbles growth rates Figs. 10 and 11 show bubbles growth and departure in microchannel. At the initial growth process, inertial force is negligible and the shape of micro bubbles is round, as shown in Fig. 11. However, inertial force increases with a bubble radius. As the inertial
force is large enough to balance the surface tension at interface, large bubbles deform and break, as shown in Fig. 10. Fig. 12 shows the simulation about bubbles growth with heating density of 70 W/cm2 and 16.16 W/cm2 during the initial process of nucleation. It can be observed that the incipient heat flux influences the bubble growth significantly. At low incipient heat flux (qincipient = 16.16 W/cm2), the surface tension controls the bubble growth. At this time Marangoni heat flow induced by temperature gradient dominates the bubble interface heat transfer. At high incipient heat flux (70 W/cm2), the surface tension hardly impact the bubbles. Meanwhile the Marangoni flow is weakened due to the little temperature gradient, and the bubble growth rates controlled by high incipient heat flux are kept high values. There were some published experimental data in similar condition in Ref. [11], and the results are also plotted in the diagram, which are consistent very well with the simulated data. Fig. 12 shows the simulation results are influenced by grid size. When bubble radius is higher than 25 lm, the difference between the results increases with calculation time. 4.4. Analysis of errors In the computation, the error in the pressure jump is investigated to estimate the accuracy of the interface tracking. The relative pressure jump error is evaluated as following:
EðDPÞ ¼
jDP DPexact j DPexact
ð26Þ
where DP is the numerical jump in pressure through bubble interface. The exact jump in pressure across a bubble is given by:
DPexact ¼ rj
ð27Þ
With two locally refined grids of h = 5 lm (Rmin/h = 2), h = 2 lm (Rmin/h = 5) in computational region, the number of total grids is 1331,817, 2593,674 respectively, and errors are plotted in Fig. 13(a). With the uniform grids of h = 15 lm (Rmin/h = 2), Fig. 10. Bubbles departure and coalescence at the second stage on the vertical symmetric profiles for q00 = 70 W/cm2. h = 20 lm, Rmin/h = 2.
Fig. 11. Axisymmetric and round bubbles (R < 40 lm) for q00 = 25 W/cm2. h = 5 lm, Rmin/h = 2.
Fig. 12. Comparing the growth of bubble radius resulted from simulation with the experimental data in Ref. [10].
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
509
Fig. 14. Volume fraction of vapor at the outlet of the channel (h = 15 lm, Rmin/h = 2).
16.16 W/cm2), bubble grows slowly, and the volume fraction of vapor at exit is low as well. 4.6. Bubble departure diameter and macro-to-micro transition threshold
Fig. 13. Relative error in calculation of pressure during bubbles growth for q00 = 70 W/cm2. (a) h = 5 lm (Rmin/h = 2), h = 2 lm (Rmin/h = 5); (b) h = 15 lm (Rmin/h = 2), and h = 20 lm (Rmin/h = 2).
h = 20 lm (Rmin/h = 2) in the computational domain, the number of total grids is 1036,298, 576,584 respectively, then errors are plotted in Fig. 13(b). Fig. 13(a) shows the errors in pressure jump, which increase with the grid size and bubble radius. Micro bubble requires much fine grids in order to obtain a solution with reasonable accuracy due to the significant refinement dependency. In addition, when a bubble grows, errors in curvature increase with time. Therefore, errors in pressure jump at interface depend primarily on grid size and time of bubble growth. Moreover, spurious velocities result from errors in curvature estimation. In Fig. 13(b), for a large bubble, errors increase quickly with bubble radius. With the bubble radius increase, drag force of liquid influences the shape of a bubble significantly, and the bubble is deformed easily, as shown in Fig. 10. Therefore, for the large deformed bubbles the refined grids may not increase the accuracy efficiently inside the wobbling region. 4.5. Vapor volume fraction at outlet Fig. 14 shows the variation of volume fraction of vapor at the outlet. The fluctuation comes from the bubble movement in the channel. There is a high value of the mean volume fraction of vapor at the outlet, it is apparent that high heat flux induces the high bubble growth rates. At low incipient heat flux (for example,
In simulation, the flow was supposed to be 1 atm at the entrance of the channel, and the velocity is 0.65 m/s. The departure diameter is 120 lm based on the criteria from Ref. [15], which is less than the channel dimension. In Figs. 6 and 10, there is no apparently jamming with one single bubble on the channel cross section. Although some small bubbles have collected into several large bubbles, their shapes are the truncated sphere still. The bubbles grow gradually due to heating in the flow. Decreasing uout is helpful to their expanding rapidly [15]. There should be the minimum value of uout for each bubble, below the value the bubble could grow beyond the confines of the micro-channel, and its shape would no longer be treated as a truncated sphere. In this simulation, the minimum uout is 0.35 m/ s, which is lower than the mean flowing velocity (0.65 m/s) at the entrance. The bubble behavior takes an important role in flow and heat transfer in micro-channel. It is proposed that the transition criterion of macro-to-micro flow might be related to the bubble departure diameter in two-phase area [1]. That is, while the departure bubble could expand to fill the channel, the further growth would be confined by the channel, the flow should be micro two-phase flow. In this paper, the dbub is less than the size of micro-channels, the shape of the single bubble is a truncated sphere, not the confined elongated bubble. 5. The simulation on the nucleate boiling in a micro trapezoid channel 5.1. Boundary conditions and source term With etching on silicon, the trapezoid channel is very popular, some work on micro-channels are on this kind of channel. For verifying the simulation with above model and method, it is shown that some calculated results are compare with the published data in Ref. [13]. The boundary conditions and fluid properties are quoted from the data in the reference. The micro-channel is 28 mm in length and the cross section is shown in Fig. 15. The
510
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
Fig. 15. The cross section of a single trapezoid micro-channel.
top of the channel is adiabatic, and the uniform heat flux is distributed at the bottom of the channel. The wall temperature can be obtained from Eq. (11). The incipient heat flux q00ONB in experiment [13] is higher than the calculated q00m , so the nucleate boiling should be controlled by the incipient heat flux q00ONB and the heat and mass rates are derived from Eqs. (21) and (23). The computational domain is arranged with locally refined grids, the grid size h is 2 lm (Rmin/h = 4) in refined domain, the number of total grids is 1092,905. 5.2. Results and discussion 5.2.1. The simulation for bubble nucleation and growth Fig. 16 displays the sequential simulated results reflecting the bubbles nucleation and their evolution at the bottom surface at the distance of 2 mm from the outlet. Here, the bubble evolution includes the process of bubble nucleation, growth and departure. Most of the bubbles are not confined from the channel dimension in this stage. Since the dimensions of the channel are different in width and height, the bubble departure diameter is properly greater than the height of the channel, that is, the bubble will be restricted in the vertical dimension in the growing process. Fig. 17 displays the bubble growth and the confined bubble on the vertical section.
The bubbles nucleate from the vicinity of the inlet to the outlet. At isotropic stage bubbles grow with a spherical shape. When the radii of bubbles reach the value of 20–25 lm, they begin to be distorted and become anisotropic due to the limitation of channel depth and the dragging of the liquid. The non-uniform growth rates among bubbles are not surprising as the local flow field may be changing due to the perturbation of the growth of other bubbles at other sites. Fig. 18(a) displays the simulated bubbles radii on position 4, 5, 6 in Fig. 16, which are compared with the experimental data in Ref. [13]. The results for bubbles on position 2, 3 while the heating density q00 is equal to 255 kW/m2 are compared with the experiments in Fig. 18(b). In the experiment, five bubble evolutions were observed on different points for different heat flux. Each process was named as cycle in the figure. The critical radius for the unconfined bubbles is approximately 22–25 lm in terms of simulation, they are closed to the values 25–30 lm from experiments [13]. Furthermore, the three simulated growth rates of bubbles (4)–(6) are the same as each other, and agree well with the referred tested data. The differences among the three positions along the flow are mainly due to the tiny differences on the enthalpy on the points. Fig. 19 shows relative errors in pressure jump. At the initial time of nucleation, bubble shape is controlled by surface tension, and errors increase slightly. When bubble boundary is close to the top wall of channel, due to limit of the channel wall and drag force of liquid, the bubble shape becomes distorted, and then errors increase quickly. It is concluded that when bubbles are disturbed by additional force, relative errors in pressure jump can not be used to estimate the accuracy of the interface tracking. Fig. 20 shows the variation of the calculated average bubbles radii with q00 = 255 kW/m2 and q00 = 189 kW/m2, which are compared with the experimental data in literature. According to the simulation, the heat flux influences the bubbles’ incipient radii and growth rates. However, it is not conspicuous in the experimental results. The experimental data in both heat fluxes correspond to the simulated profile with 255 kW/m2 of heat flux in the growing period. Meanwhile, in the initial nucleating period, the observed bubble radii are larger than the simulated values.
Fig. 16. Bubbles nucleation and growth at the bottom surface downstream for one cycle (t = 0–140 ms), G = 341 kg/m2, q00 = 189 kW/m2, h = 2 lm, Rmin/h = 4.
Fig. 17. Bubbles nucleation and growth at the vertical wall upstream (t = 0–300 ms), G = 341 kg/m2 s, q00 = 189 kW/m2, h = 2 lm, Rmin/h = 4.
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
511
Fig. 20. Time evolution of average bubble radius (h = 2 lm, Rmin/h = 4).
Fig. 18. Time evolution of bubble radius for G = 341 kg/m2 s: (a) q00 = 189 kW/m2, (b) q00 = 255 kW/m2; h = 2 lm, Rmin/h = 4. Fig. 21. Pressure variation at the horizontal section near the outlet (h = 2 lm, Rmin/h = 4).
the maximum pressure that can be sustained in the bubble is dictated significantly by the saturation pressure corresponding to the wall temperature. As shown in Fig. 21, the pressure spike occurs at the nucleation sites which cause the pressure variation along the flow direction. The pressure jump at the nucleation site leads to a disruption in the flow. When bubbles is sticking on the wall, the surface tension force along the contact line may balance the drag force of the incoming liquid, and can disturb flow at interface. 6. Summary
Fig. 19. Relative error in calculation of pressure during bubbles growth for G = 341 kg/m2 s, q00 = 189 kW/m2, h = 2 lm, Rmin/h = 4.
The nucleate boiling in micro-channels was investigated with numerical simulation, and some simulation results were compared with previous experiments in literature. The comparison indicates that the simulation is feasible and reliable. The conclusions are as following:
5.2.2. Flow stability during nucleate boiling in a micro-channel According to Eq. (10), the pressure jump occurs at a bubble interface due to the surface tension. As a bubble grows on a cavity,
(1) Through the analysis with Marangoni heat transfer on bubble surface, it is proposed that Marangoni heat transfer is important at the onset of nucleate boiling in micro-channels. Moreover, the nucleate boiling process is divided into two
512
R. Zhuan, W. Wang / International Journal of Heat and Mass Transfer 53 (2010) 502–512
stages. In the first stage, bubbles are suppressed by surface tension, and Marangoni convection dominates the interface heat transfer. In the second stage, the wall superheat temperature is higher than the transitional point, and qONB is large enough to overcome the surface tension. At this time, bubbles grow rapidly and lift off from the wall. (2) The feature of nucleate boiling in common channel is free bubble flow, while in micro-channels it is a confined bubble flow. It is proposed that the macro-to-micro transition criterion in two-phase flow and heat transfer should be related to the bubble departure diameter, furthermore, the width and height might influence the transitional point too. The departure diameter could be considered as the transitional criterion of macro-to-micro dimension. (3) It is proved that the incipient heat flux would affect the bubble incipient radius and the bubble growth rates. In addition, the higher the heat flux, the higher the volume fraction of vapor at the outlet of a micro-channel. (4) The pressure spike occurs at the nucleation sites due to the influence of surface tension on interface. The pressure perturbation at the nucleation sites might be the causes of disruption in the flow.
Acknowledgement The authors appreciate the support from National Natural Science Fund in China (50576054). References [1] John R. Thome, Boiling in microchannels: a review of experiment and theory, Int. J. Heat Fluid Flow 25 (2004) 128–139. [2] S.G. Kandlikar, Fundamental issues related to flow boiling in minichannels and microchannels, Exp. Thermal Fluid Sci. 26 (2002) 389–407. [3] A.E. Bergles, V.J.H. Lienhard, G.E. Kendall, P. Griffith, Boiling and evaporation in small diameter channels, Heat Transfer Eng. 24 (2003) 18–40. [4] S.M. Ghiaasiaan, R.C. Chedester, Boiling incipience in microchannels, Inter. J. Heat Mass Transfer 45 (2002) 4599–4606. [5] R.C. Chedester, S.M. Ghiaasiaan, A proposed mechanism for hydrodynamicallycontrolled onset of significant void in microtubes, Int. J. Heat Fluid Flow 23 (2002) 769–775. [6] J.E. Kennedy, G.M. Roach Jr., M.E. Dowling, S.I. Abdel-Khalik, S.M. Ghiaasiaan, S.M. Jeter, Z.H. Qureshi, The onset of flow instability in uniformly heated horizontal microchannels, ASME J. Heat Transfer 122 (2000) 118–125. [7] Y.Y. Hsu, On the size range of active nucleation cavities on a heating surface, J. Heat Transfer 84 (1962) 207–216. [8] S.G. Kandlikar, V.R. Mizo, M.D. Cartwright, Ikenze, Bubble nucleation and growth characteristics in subcooled flow boiling of water, HTD-vol. 342, in: ASME Proceedings of the 32nd National Heat Transfer Conference, vol. 4, 1997, pp. 11–18. [9] W. Qu, I. Mudawar, Prediction and measurement of incipient boiling heat flux in microchannel heat sinks, Inter. J. Heat Mass Transfer 45 (2002) 3933–3945. [10] J. Li, P. Cheng, Bubble cavitation in a microchannel, Int. J. Heat Mass Transfer 47 (2004) 2689–2698. [11] Dong Liu, Poh-Seng Lee, Suresh V. Garimella, Prediction of the onset of nucleate boiling in microchannel flow, Int. J. Heat Mass Transfer 48 (2005) 5134–5149. [12] N.O. Young, J.S. Goldstein, M.J. Block, The motion of bubbles in a vertical temperature gradient, J. Fluid Mech. 6 (1959) 350–356. [13] P.C. Lee, F.G. Tseng, Chin Pan, Bubble dynamics in microchannels. Part I: single microchannel, Int. J. Heat Mass Transfer 47 (2004) 5575–5589. [14] H.Y. Li, F.G. Tseng, Chin Pan, Bubble dynamics in microchannels. Part II: two parallel microchannels, Int. J. Heat Mass Transfer 47 (2004) 5591–5601. [15] Weilin Qu, Issam Mudawar, Prediction and measurement of incipient boiling heat flux in micro-channel heat sinks, Int. J. Heat Mass Transfer 45 (2002) 3933–3945.
[16] S.S. Mehendal, A.M. Jacobi, R.K. Shah, Fluid flow and heat transfer at micro- and meso-scales with application to heat exchanger design, Appl. Mech. Rev. 53 (7) (2000) 175–193. [17] P. Kew, K. Cornwell, Correlations for prediction of boiling heat transfer in small diameter channels, Appl. Thermal Eng. 17 (8–10) (1997) 705– 715. [18] Sanja Petrovic, Tony Robinson, Ross L. Judd, Marangoni heat transfer in subcooled nucleate pool boiling, Int. J. Heat Mass Transfer 47 (2004) 5115– 5128. [19] J. Li, G.P. Peterson, Microscale heterogeneous boiling on smooth surfaces— from bubble nucleation to bubble dynamics, Int. J. Heat Mass Transfer 48 (2005) 4316–4332. [20] J.L. McGrew, F.L. Bamford, T.R. Rehm, Marangoni flow: an additional mechanism in boiling heat transfer, Science 153 (1966) 1106–1107. [21] P. Arlabosse, L. Tadrist, H. Tadrist, J. Pantaloni, Experimental analysis of the heat transfer induced by thermocapillary convection around a bubble, J. Heat Transfer 122 (2000) 66–73. [22] W.F. Noh, P.R. Woodward, SLIC (Simple Line Interface Method), in: Lecture Notes in Physics, vol. 59, Springer, 1976. [23] C.W. Hirt, B.D. Nichols, Volume of Fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 39 (1981) 201–225. [24] D.L. Youngs, Time-dependent multi-material flow with large fluid distribution, in: K.W. Morton, M.L. Norman (Eds.), Numerical Methods for Fluid Dynamics, 1986, pp. 187–221. [25] N. Ashgriz, J.Y. Poo, FLAIR: flux line-segment model for advection and interface reconstruction, J. Comput. Phys. 93 (1991) 449–468. [26] S.O. Kim, H.C. No, Second order model for free surface convection and interface reconstruction, Int. J. Numer. Methods Fluids 26 (1) (1998) 79–100. [27] J.E. Polliod, An analysis of piecewise linear interface reconstruction algorithm for volume of fluid methods, M.Sc. Thesis, Department of Mathematics, University of California, Davis, 1992. [28] W.J. Rider, D.B. Kothe, Reconstructing volume tracking, J. Comput. Phys. 141 (1998) 112–152. [29] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (1992) 335–354. [30] J.U. Brackbill, D.B. Kothe, Dynamic modeling of the surface tension, in: Proceedings of the 3rd Microgravity Fluid Physics Conference, Cleveland, OH, 1996, pp. 693–698. [31] M. Seifollahi, E. Shirani, N. Ashgrizb, An improved method for calculation of interface pressure force in PLIC-VOF methods, Eur. J. Mech. B/Fluids 27 (2008) 1–23. [32] E. Shirani, N. Ashgriz, J. Mostaghimi, Interface pressure calculation based on conservation of momentum for front capturing methods, J. Comput. Phys. 203 (2005) 154–175. [33] Daniel Lörstad, Laszlo Fuchs, High-order surface tension VOF-model for 3D bubble flows with high density ratio, J. Comput. Phys. 200 (2004) 153– 176. [34] Marianne M. Francois, Sharen J. Cummins, Edward D. Dendy, Douglas B. Kothe, James M. Sicilian, Matthew W. Williams, A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework, J. Comput. Phys. 213 (2006) 141–173. [35] H. Wang, X.F. Peng, W.K. Lin, C. Pan, B.X. Wang, Bubble-top jet flow on micro wires, Int. J. Heat Mass Transfer 47 (2004) 2891–2900. [36] H. Wang, X.F. Peng, B.X. Wang, D.J. Lee, Bubble-sweeping mechanisms, Sci. China, Ser. E 46 (2003) 225–233. [37] Hao Wang, Xiaofeng Peng, Suresh V. Garimella, David M. Christopher, Microbubble return phenomena during subcooled boiling on small wires, Int. J. Heat Mass Transfer 50 (2007) 163–172. [38] J.R. Thome, V. Dupont, A.M. Jacobi, Heat transfer model for evaporation in microchannels. Part I: presentation of the model, Int. J. Heat Mass Transfer 47 (2004) 3375–3385. [39] Nengli Zhang, David F. Chat, Models for enhanced boiling heat transfer by unusual Marangoni effects under microgravity conditions, Int. Comm. Heat Mass Transfer 26 (1999) 1081–1090. [40] Satish G. Kandlikar, Nucleation characteristics and stability considerations during flow boiling in microchannels, Exp. Therm. Fluid Sci. 30 (2006) 411– 417. [41] C.L. Vandervort, A.E. Bergles, M.K. Jensen, Heat transfer mechanisms in very high heat flux subcooled boiling, ASME Fundam. Subcooled Flow Boiling HTD 217 (1992) 1–9. [42] E.J. Davis, G.H. Anderson, The incipience of nucleate boiling in forced convection flow, AIChE J. 12 (1966) 774–780. [43] Jr. Hung Tsai, Liwei Lin, Transient thermal bubble formation on polysilicon micro-resisters, J. Heat Transfer 124 (2002) 375–382. [44] Liwei Lin, Microscale thermal bubble formation: thermophysical phenomena and applications, Microscale Thermophys. Eng. 2 (1998) 71–85.