Strojniški vestnik - Journal of Mechanical Engineering 51(2005)7-8, 436-444 UDK-UDC 536.2 Izvirni znanstveni članek - Original scientific paper (1.01) Simulation of Void Fraction Profile Evolution in Subcooled Nucleate Boiling in a Vertical Annulus with a Bubble-Tracking Model Ivo Kljenak*, Borut Mavko Reactor Engineering Division, Jozef Stefan Institute, Jamova 39, Ljubljana, Slovenia,* ivo.kljenak@ijs.si Abstract A three-dimensional bubble-tracking model of subcooled nucleate boiling flow in a vertical channel at low-pressure conditions is proposed, with specific application to the case of boiling in an annulus with a central heating rod. In the model, vapour is distributed in the liquid in the form of individually tracked bubbles. The overall behaviour of the liquid-vapour system results from motion, interaction, coalescence and boiling mechanisms prescribed mostly at the level of bubbles. The wall heat transfer coefficient and the wall temperature are calculated from one-dimensional correlations. The partitioning of the heat flux, which is consumed for bubble nucleation and heating of the liquid, varies along the flow and depends on bubble size as well as on local flow conditions. Bubbles are nucleated with constant frequencies at fixed nucleation sites randomly distributed over the heated surface. Liquid temperature profiles at different axial locations are determined from steady-state energy balances. The nucleation site density is determined from a balance between vapour generation rate, bubble departure sizes and nucleation frequencies. After nucleation, bubbles slide on the heated surface, detach and then gradually migrate into the low-temperature region away from the heated surface, where they eventually condense. Both bubble detachment and migration are modelled probabilistically. Bubble lateral migration is restricted by the lift force due to the liquid velocity gradient. The proposed model was applied to experiments on subcooled boiling that were carried out at Purdue University (USA) by Bartel [1]. A good agreement between measured and calculated void fraction profiles at different axial locations was obtained. Introduction In the subcooled part of upward nucleate boiling flow in a vertical channel with a heated wall, the temperature near the wall and the bulk fluid temperature are respectively higher and lower than the saturation temperature. Subcooled boiling is thus characterized by a "higher-temperature" two-phase region near the heated surface and a "lower-temperature" single-phase liquid region away from the heated surface. The evolution of void fraction in subcooled boiling flow may be modelled using various approaches with different time and length scales. One-dimensional two-fluid models with various degrees of empiricism may predict fairly well the void fraction, averaged over the channel cross-section (Hari and Hassan [2], Končar and Mavko [3]). On the other hand, models based on local instantaneous description of the flow are at present still computationally too demanding to be applied to boiling systems which may have a complex interface structure due to the presence of up to several thousand bubbles. "Intermediate-level" models, which may be applied to nucleate boiling flows, include multidimensional two-fluid models, based on ensemble and volume averaging of local instant conservation equations (Kurul and Podowski [4], Janssens-Maenhout et al. [5], Končar et al. [6], Lee et al. [7]), and so-called bubble-tracking models, in which gas is distributed in the liquid in the form of individually-tracked bubbles (Mortensen and Trapp [8]). Among experimental results on subcooled boiling in channels, not many authors have measured the non-homogeneous radial distributions of two-phase flow parameters, such as void fraction and bubble size. Recently, these kinds of experiments have been carried out by Bartel [1] and Lee et al. [7]. In the present work, a three-dimensional bubble-tracking model of subcooled nucleate boiling flow in a vertical channel is presented. The behaviour of the bubble 436 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 Nomenclature Greek letters A channel cross-section [m2] a void fraction [-] Dh channel hydraulic diameter [m] p density [kg/m3] G mass flux [kg/m2s] G surface tension [N/m] Ro radius of annulus outer wall [m] standard deviation of bubble diam. dist. [m] Ri radius of annulus inner wall [m] U thermal diffusivity [m2/s] T temperature [K] H viscosity [Pa·s] V kinematic viscosity [m2/s] c specific heat [J/kg-K] d bubble equivalent diameter [m] f bubble nucleation frequency [s-1] Subscripts g gravitational acceleration [m/s2] b bubble h heat transfer coefficient [W/m2K] bl bubble relative velocity specific enthalpy [J/kg] d bubble departure equivalent diameter j volumetric flux [m/s] g gas k thermal conductivity [W/m-K] i annulus inner wall l bubble vertical chord length [m] i-th bubble p pressure [Pa] l liquid probability [-] loo undisturbed liquid velocity q'' heat flux [W/m2] o annulus outer wall r distance from annulus inner wall [m] p constant pressure t time [s] sat saturation conditions w velocity in z-direction [m/s] w heated wall z axial coordinate [m] Prl liquid Prandtl number = vl / v>l Other symbols Reb bubble Reynols number = wbl · db /vl < > average over channel cross-section population is simulated by considering each bubble separately. The overall behaviour of the liquid-vapour system results from motion, interaction and boiling mechanisms prescribed mostly at the bubble level. The proposed work represents a further development of a model, which has already been presented earlier. The model was first developed for subcooled boiling in a cylindrical tube (Kljenak [9]). It was later extended to annular channels, and calculated results were compared to experimental measurements at a single axial location (Kljenak et al. [10], Kljenak et al. [11]) which were obtained at Seoul National University (Lee et al. [7]). In the present work, the model was used to simulate experiments with boiling water in a heated vertical annular channel at atmospheric pressure, which were performed by Bartel [1]. In these experiments, flow parameters at different axial locations along the flow were measured. Thus, experimental and simulated evolutions of radial void fraction profiles are compared. Physical model In the present work, the equations apply to boiling flow in an annular channel. However, the model may be applied to a cylindrical tube as well. Bubble axial motion and interaction In the proposed model, bubbles assume a rigid ellipsoidal shape and move upwards with their symmetry axis always vertical. The velocity of a bubble is calculated by first adding the bubble relative velocity (calculated from a correlation by Peebles and Garber, as cited by Wallis [12]) to the local hypothetical “undisturbed” liquid velocity, obtained from the 1/7th power law. Then, if some other (leading) nearby bubble is found to be present ahead of the bubble whose velocity is being calculated, an increase due to wake drift is added. The liquid velocity behind bubble i, which is increased due to wake drift, is described in, basically, the same form as suggested by Bilicki and Kestin [13]: wl (z, r, t) = wlarj (z, r, t) + 2g/3 lb max 2(zi-z) (wbi (t) - wlao (zi, ri, t)) (1) where lb max denotes the bubble maximum vertical chord length (ellipsoid vertical axis) and t, is an empirical attenuation factor, set equal to 3.0, which was introduced to obtain a weaker wake drift as the velocity in the wake decreases with distance from the wake axis. Simulation of void fraction profile evolution in subcooled nucleate boiling 437 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 A necessary condition for a bubble to be influenced by a leading bubble through wake drift and eventually collide with it is that bubbles overlap laterally (Fig. 1) by more than a certain critical fraction, called minimum relative overlapping. A similar approach was already proposed by Mortensen and Trapp [8]. If, following axial collision, bubbles overlap in space, the upper bubble is displaced laterally (sideways) for a fraction of its width if there is no other nearby bubble to prevent the movement. Otherwise, the upper bubble is displaced upwards so that bubbles barely stick. If bubbles still stick after collision (that is, if they were not separated due to a lengthy lateral displacement of the upper bubble), they remain sticking, move along together with the upper bubble's velocity and eventually merge if they do not separate earlier due to either turbulent dispersion, subsequent movements of either bubble or coalescence of either bubble with some other bubble. If separate bubbles do not overlap laterally more than the critical fraction, the motion of the trailing bubble is not affected by the leading bubble and bubbles behave as if they would not overlap at all. This rule was prescribed to approximate the influence of bubble agitation, which occurs in real bubbly flow and allows tightly packed bubbles to overtake one another. The drawback of this approach is that bubbles may briefly overlap in space, which is not physically realistic. lower bubble upper bubble Prodanovic et al. [16] observed that bubbles usually slide a couple of diameters before detaching. In the proposed model, the following empirical approach was adopted, based on above-cited experiments: after nucleation, bubbles slide along the heated surface for some distance before attempting detachment. During sliding, the ellipsoidal bubble vertical axis is longer than the bubble horizontal axis to approximate the bubble inclination which was observed in experiments. Detachment, which may occur if no other bubble obstructs the radial motion away from the heated surface, is modelled probabilistically in the same way as bubble radial migration (see below). After detaching from the heated surface, bubbles tend to migrate towards the lower-temperature region (again, if other bubbles do not obstruct their motion), where they eventually condense. After detachment, the bubble shape changes and the horizontal axis is longer than the vertical axis. In bubbly flow in general (that is, in boiling as well as in adiabatic flow), bubbles are distributed over the channel cross-section, supposedly as the result of interaction of different phenomena: liquid turbulent flow, transverse lift force, wall lubrication force and bubble interaction (Žun [17], Liu [18], Ohnuki and Akimoto [19], Okawa et al. [20]). Bubble transverse motion over the channel cross-section is presumably partly random, due to the interaction of bubbles with turbulent eddies. A probabilistic approach was thus implemented to model bubble radial motion (migration) towards the tube outer wall: the motion consists of finite steps that are equal to a fraction of bubble width, each displacement occurring with a certain probability. At present, the proposed model was developed for boiling systems in which all bubbles are located between the heated inner annulus wall and the middle of the annular gap. The lift force, which is related to the liquid velocity gradient, is assumed to represent a restraining force to bubble radial motion away from the heated wall. Thus, the probability of migration pm was modelled to increase with decreasing velocity gradient over the channel cross-section: Figure 1. Relative overlapping between lower and upper bubble = s/b The two-dimensional “undisturbed” liquid velocity profile causes bubbles located at different distances from the wall to move with different velocities, thus promoting bubble collisions and coalescence. Bubble sliding, detachment and radial motion In various experiments on subcooled nucleate boiling of water at low-pressure conditions, it has been observed that bubbles nucleated on a heated wall first slide along the wall and then tend to detach and migrate towards the tube region away from the heated surface (Bibeau and Salcudean [14], Zeitoun and Shoukri [15], Prodanovic et al. [16]). To the authors’ knowledge, no consensus concerning bubble sliding distance has been reached yet. pm =1.0-Cr ?wl l? ?r 1/2 (2) where Cr is an empirical coefficient. As bubble lateral motion may be affected by turbulent eddies of a comparable size as the bubble, bubble migration is attempted every time a bubble moves axially a distance equal to its maximum vertical chord length. Turbulent dispersion The relative motion between bubbles is mainly influenced by eddy motion of the length scale of bubble size (Prince and Blanch [21]). The influence of turbulent eddies, which may affect wake drift or sticking bubbles, is modelled as a succession of random binary events. Thus, each event may have two possible outcomes: at a 438 Kljenak I. - Mavko B. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 given instant, bubble motion is or is not disturbed by turbulent dispersion. Higher turbulence intensity is simulated by prescribing a higher probability of dispersion. Intensities of turbulent dispersion are assumed to be constant over the channel cross-section. Disruption of wake drift and of sticking bubbles by turbulent eddies is related to the turbulence length scale and is simulated every time a bubble has moved in the axial direction a distance equal to 1/20 of the channel hydraulic diameter. Bubble coalescence Following axial collision and sideways or upwards displacement of the upper bubble, bubbles which still overlap more than the minimum relative overlapping stick together. Bubbles eventually merge after sticking together for a certain time interval (so-called "rest time") if they are not dispersed by liquid turbulence or do not move apart due to axial collisions with other bubbles. Bubble coalescence occurs instantly only if the leading bubble is sliding on the heated surface whereas the trailing bubble is not, as the impact between bubbles is presumably stronger due to larger velocity differences. In the present work, it was assumed that the impact following bubble lateral collision is not strong enough to cause rupture of the vapour-liquid interface, so that coalescence following bubble lateral collision was not modelled. from thermal energy balances. The gas phase is assumed to be at saturation conditions. The temperature of the heated wall is assumed to increase until it reaches a value determined from a correlation by Shah (as cited by Kandlikar [23]): q'' = (230(Gh lg ) 0.5hw(Tw -Tsat))2 (5) where the single-phase heat transfer coefficient h1t, has to be calculated from the well-known Dittus-Boelter correlation (Collier [24]) and hlg indicates the difference between vapour and liquid specific enthalpies at saturation conditions. Before reaching that value, the wall temperature is calculated from the relation: q''= h10(Tw- < Tl >) (6) where denotes the liquid temperature, averaged over the channel cross-section. Partitioning of wall heat flux In the proposed model, the wall heat flux is partitioned as follows: q'' = Cl h1^(Tw - ) + q''nucl + q'' ' slid (7) Liquid temperature In the proposed model, the liquid temperature Tl depends on the distance from the inner heated wall r and obeys the following law (Sekoguchi et al. [22]): Ti -Tl Ti - To R o - Ri 1/m (3) where the exponent m may depend on the flow rate. Namely, at higher flow rates, the temperature gradient near the wall is expected to be somewhat steeper due to more intense turbulent mixing. Liquid temperature profiles at different axial locations along the tube are obtained using steady-state values of the average cross-sectional enthalpy . The liquid temperature profile must be such that the liquid specific enthalpy hl fulfils the condition: (h\ \((1 - ?)? l + ?? g)dA = A \(1-?)?lhldA + ??ghgdA AA (4) where ? denotes the local void fraction and integrals are calculated over the channel cross-section. Steady-state values of at different axial locations are obtained The first term on the r.h.s. of Eq. (7) represents heat transfer due to single-phase forced convection. The factor Cl accounts for the portion of the heated surface not covered by bubbles. The term q''nucl denotes the heat flux consumed for bubble nucleation. The term q''slid, which denotes the heat flux consumed for growing of bubbles that slide on the heated surface, is determined as in the work of Tsung-Chang and Bankoff [25]: q'' slid = 2kl(Tw -Tsat) (8) The surface through which heat is transferred to the bubble is represented by a circular area with a diameter equal to the bubble vertical axis. Bubble nucleation In the proposed model, bubbles are nucleated at fixed nucleation sites randomly distributed over the heated surface, instantly reach departure size and assume an ellipsoidal shape. The bubble equivalent departure diameter dd is constant at each site. Bubble diameters over nucleation sites are distributed according to Gaussian distributions and are randomly generated in intervals [dd - 3a, dd + 3a], where dd denotes the local mean bubble equivalent diameter and o the local standard deviation. The mean bubble size was assumed to depend on local subcooling and was calculated in the same way r Simulation of void fraction profile evolution in subcooled nucleate boiling 439 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 as in the work of Končar et al. [6], using Unal’s [25] mechanistic model. Unal’s model describes the bubble departure diameter dd as a function of pressure, liquid subcooling, heat flux and liquid flow velocity: 2.42-10 5 ¦ p0709a d = Cbw---------------b=------------- ¦\]txp dd=C (9) where coefficients a, b and are defined as: a = (q''-h 1^ ¦ ATsub)13k 2C1 3hlg Pg^l™ljPlcpl V klPl lpl C = b = h lg^'[cpl/(0.013hlgPrl1.7)]3 Ja/(Pl-pg)g ATsub 2(1-PgI Pl) = vl 0.61 0.47 for vl > 0.61 m/s for vl < 0.61 m/s (10) (11) (12) (13) The range of applicability of the correlation is: - pressure: 0.1 < p < 17.7 MPa, - wall heat flux: 0.47 < q'' < 10.64 MW/m2, - liquid velocity: 0.08 < wl < 9.15 m/s, - liquid subcooling: 3.0 < ATsub < 86 K. Since, in the present work, the heat flux in the considered experimental data is below the range of applicability of the correlation, the coefficient Cbw was added in Eq. (9) to describe relatively large bubbles at low-pressure conditions. The frequency of bubble nucleation at individual sites is calculated from a correlation by Cole (1960, as cited by Ivey [27]): f = 4g(J>i-pg) 3dd pt 1/ 2 (14) Bubble evaporation and condensation In the proposed model, bubbles may further grow while part of them is still within the region near the heated wall where the temperature is higher than the saturation temperature. The interfacial heat transfer coefficient hint is calculated from a correlation already used by Mortensen and Trapp [8]: db (2 (0.4Re 0.06Re b 0.67 )Prl 0.4 (15) Bubbles that move laterally into the lower-temperature region collapse instantly if the liquid temperature at the bubble tip closest to the heated surface is lower than the saturation temperature. Numerical model Bubble behaviour The proposed model was implemented as a computer code. Bubble axial motion is simulated with a simple discrete time-step method, neglecting inertial effects: zi (t + At) = zi (t) + wbi (t) ¦ At (16) As bubbles in the proposed model undergo significant accelerations only briefly before axial collision with a leading bubble or after radial migration to a higher liquidvelocity region, the added mass effect is not taken into account. After each axial displacement during a time step, bubbles, which overlap more than the minimum relative overlapping, are adjusted if they also overlap in space. Adjustments start at the tube entrance, and upper bubbles are adjusted with respect to lower bubbles. If possible, each upper bubble moves laterally (sideways) for up to a fraction of its width. If this is not possible due to the presence of other bubbles, the upper bubble is displaced upwards (see Section “Bubble axial motion and interaction”). This adjustment of bubbles simulates bubble collisions and subsequent displacements of upper bubbles. Bubbles' cross-sectional coordinates assume discrete values, which correspond to points located on concentric circles, centred on the annulus axis (Fig. 2). The distance between neighbouring points along concentric circles is constant. The distance between any neighbouring points must be of the order of a fraction of the smallest bubbles' width (usually about few hundredths). Bubble lateral movements are modelled as instantaneous jumps to other points and occur between time steps. Mergers between bubbles, bubble condensations and bubble nucleations are also modelled as instantaneous events which occur between time steps. 1 440 Kljenak I. - Mavko B. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 Figure 2. Schematic representation of bubble centres' discrete cross-sectional coordinates. Volume fraction and energy coupling The vertical channel is divided in the axial direction into control volumes. The partitioning of the heat flux into vapour generation and liquid heating is constant within each control volume. The total simulation time is divided into time sub-intervals during which passages of bubbles through control volume boundaries are recorded. At the end of each sub-interval, the liquid velocity profile corresponding to each boundary is adjusted to satisfy the following mass conservation condition (Kowe et al. [28]): Ajl= \wl? (1-?)dA +Cmwbl ?dA AA ) where ? denotes the time-averaged local void fraction obtained from recordings of bubble passages at the control volume boundary and integrals are calculated over the channel cross-section. The coefficient of added mass Cm was set equal to 0.5. In the same way, the liquid temperature profile is set so that the cross-sectional enthalpy at each boundary assumes the steady-state value (Eq. 4), obtained from the total heat input below the boundary. These liquid velocity and temperature profiles are then used in calculations during the next time sub-interval. Within each control volume, local liquid velocities and temperatures are calculated by linear interpolation between values, which correspond to the control volume lower and upper boundaries. These values are needed when calculating liquid velocities and temperatures corresponding to bubble centres. For each control volume, the partitioning of the heat flux into vapour generation and liquid heating is based on conditions at the control volume upper boundary and is determined periodically after each time sub-interval. The nucleation site density at a given axial location is determined by an iterative calculation so that the sum of all the heat flux components equals the imposed heat flux. Nucleation sites on the tube wall are placed randomly and used to generate bubbles during the next time sub-interval. Results Experimental conditions Experiments on subcooled nucleate boiling at atmospheric pressure were performed at Purdue University (USA) by Bartel [1]. The experiments were carried out in a vertical annulus with a heated inner rod. The diameter of the rod was 19.1 mm, whereas the inner diameter of the outer tube was 38.1 mm. The length of the heated part of the annulus was 1.5 m. Local void fraction was measured using an electrical conductivity probe technique. The data were collected simultaneously at different axial locations. Thus, the axial evolution of radial distribution of flow parameters was observed. Experimental conditions for runs, which were simulated with the proposed model, are presented in Table 1. The inlet subcooling refers to conditions at atmospheric pressure. Table 1. Experimental conditions Run ?Tsub [°C] q’’ 2 [kW/m2] G [kg/m2s] 1 8.9 105 470 2 6.1 128 701 3 4.8 128 701 4 5.2 145 700 Model parameters At each nucleation site, the bubble departure diameter was constant but generated as a random variable with a Gaussian distribution (see section “Bubble nucleation”). It was assumed that the smallest bubble diameter (lower boundary of the interval: dd - 3a) was always 0.0005 m, from which the interval upper boundary was calculated. The value of the coefficient Cbw in Eq. (9) was set to 1.5. The ratio of horizontal to vertical bubble axis before bubble detachment from the wall was set to 0.8, to approximate the inclination of bubbles while they are sliding on the walls. After detachment, the ratio was set to 1.2. All simulations were carried out with identical values of bubble minimum relative overlapping (0.3) and rest time during which bubbles stick together before merging (0.02 s). There is a lack of information on rest times in turbulent flows, necessitating the use of what is in effect an adjustable parameter (Prince and Blanch [21]). Bubble radial motion (including detachment from the heated wall and sideways displacement of the upper bubble after axial collision) consisted of finite steps of 1/20 bubble width. The factor Cr in Eq. (2) was set equal to 0.03 for all runs. The probabilities of turbulent dispersion were set to 0.09 for G=470 kg/m2s and 0.13 for G=700-701 kg/m2s. Simulation of void fraction profile evolution in subcooled nucleate boiling 441 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 The ratios of probabilities of turbulent dispersion correspond to the ratios of bulk Reynolds numbers. The factor m in Eq. (3) (liquid temperature profile) was set to 4.0 for all runs. For G=470 kg/m2?s, it was assumed that bubbles slide along the wall for a distance equal to twice their maximum vertical chord length (that is, twice the ellipsoid vertical axis) before attempting detachment. For G=700-701 kg/m2?s, it was assumed that the “sliding distance” is equal to one vertical axis. The rationale for the difference is that bubble detachment is more likely to occur sooner at higher mass flow rates due to the influence of turbulent eddies. 0.15 0.10 0.05 0.00 0.0 0.2 0.4 0.6 0.8 1.0 r / (Ro - Ri ) Simulation results Figures 3-6 show experimental and simulated time-averaged void fraction radial profiles at different experimental conditions and different axial locations along the channel. The coordinate z=0 corresponds to the beginning of the heated section. In general, the overall agreement between simulations and experiments is good. As bubbles instantly assume the shape of ellipsoids, which have their axis always vertical and touch the wall only with their tip, void fraction assumes zero values at the heated wall. Due to the assumption of rigid ellipsoidal bubbles, whose axis always remain vertical, the proposed approach is necessarily limited to relatively low void fraction values. Namely, at higher void fractions, bubbles in actual flows are packed more closely, so that their shape is probably somewhat distorted and their axis do not always remain vertical. In the proposed model, parameters that were adjusted either assume a constant value or their variation may be justified on physical grounds. Thus, the presented simulations of experiments, performed with different combinations of boundary conditions (subcooling, mass flux and heat flux), augur that the model should be applicable to subcooled boiling in vertical annular channels of similar dimensions and over a range of experimental conditions. Figure 3. Simulated and experimental void fraction profiles (G=470 kg/m2s, q’’=105 kW/m2,?Tsub =8.9° C) 0.20 0.15 0.10 0.05 0.00 0.0 0.2 0.4 0.6 0.8 1.0 r / (Ro - Ri ) Figure 4. Simulated and experimental void fraction profiles (G=701 kg/m2s, q’’=128 kW/m2, ?Tsub =6.1° C) 0.25 0.20 0.15 0.10 0.05 0.00 ¦r'V z/Dh exp. theory 6.9 • ............... 44.6 ? -------- 65.6 A -------- R/- \ - 7' v> D \ -»'"¦•¦ Vv \ / * \n \ ^ t .............."¦¦-- Vjj A 0.0 0.2 0.6 0.4 r / (Ro - Ri ) 0.8 1.0 Figure 5. Simulated and experimental void fraction profiles (G=701 kg/m2s, q’’=128 kW/m2, ?Tsub =4.8° C) 442 Kljenak I. - Mavko B. Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 0.1 0.0 0.0 0.2 0.4 0.6 r / (Ro - Ri ) 0.8 1.0 Figure 6. Simulated and experimental void fraction profiles (G=700 kg/m2s, q’’=145 kW/m2,?Tsub =5.2° C) Conclusions A three-dimensional bubble-tracking model, in which empiricism is included at a "more fundamental" level, was developed to simulate upward subcooled nucleate boiling flow in a vertical annular tube with a central heating rod. Bubble collective behaviour results from motion, interaction and heat transfer mechanisms prescribed at the level of individual bubbles. The model was used to simulate experiments performed with water at near atmospheric pressure. A good overall agreement between calculated and measured radial profiles of void fraction at different axial locations along the channel was obtained. Although the model contains a number of adjustable parameters, the comparison of simulated results with experimental data indicates that the proposed approach captures the basic mechanisms that govern the development and evolution of subcooled nucleate boiling along a heated channel at near atmospheric pressure and relatively low void fraction. References [1] M.D. Bartel, Experimental Investigation of Subcooled Boiling, M.Sc.Thesis, Purdue Univ., West Lafayette, IN, USA, 1999. [2] S. Hari, Y.A. Hassan, Improvement of the subcooled boiling model for low-pressure conditions in thermal-hydraulic codes, Nuclear Engng. Design 216 (2002) 139-152. [3] B. Končar, B. Mavko, Modelling of low-pressure subcooled flow boiling using the RELAP5 code, Nuclear Engng. Design 220 (2003) 255-273. [4] N. Kurul, M.Z. Podowski, On the modeling of multidimensional effects in boiling channels, Proc. 27th National Heat Transfer Conf., Minneapolis, USA, 1991. [5] G. Janssens-Maenhout, JU. Knebel, U. Mueller, Subcooled nucleate boiling at low pressure and low heat flux, Proc. 3rd Int.Conf. Multiphase Flow ICMF'98, Lyon, France, 1998. [6] B. Končar, I. Kljenak, B. Mavko, Modelling of local two-phase flow parameters in upward subcooled flow boiling at low pressure, Int.J. Heat Mass Transfer 47 (2004) 1499-1513. [7] T.H. Lee, G.C. Park, D.J. Lee, Local flow characteristics of subcooled boiling flow of water in a vertical concentric annulus, Int. J. Multiphase Flow 28 (2002) 1351-1368. [8] GA. Mortensen, J.A. Trapp, Two-phase flow modeling with discrete particles, Proc. ASME Heat Transfer Div., Two-Phase Flow in Energy Exchange Systems, vol. 220, 1992, pp. 73-85. [9] I. Kljenak, Modeling of void fraction and liquid temperature profiles evolution in vertical subcooled nucleate boiling flow, Proc. ASME-ZSITS IntThermal Science Seminar, Bled, Slovenia, 2000. [10] I. Kljenak, G.C. Park, B. Mavko, T. Lee, Bubble-tracking modeling of subcooled nucleate boiling in a vertical annulus, Proc. 4th Int. Conf. Multiphase Flow, New Orleans, USA, 2001. [11] I. Kljenak, G.C. Park, B. Mavko, T. Lee, Bubble-tracking modeling of subcooled nucleate boiling in a vertical annulus, Proc. 12th Int. Heat Transfer Conf., Grenoble, France, 2002. [12] GB. Wallis, One-Dimensional Two-Phase Flow, McGraw-Hill, 1969, pp. 248-251. [13] Z. Bilicki, J. Kestin, Transition criteria for two-phase flow patterns in vertical upward flow, Int.J. Multiphase Flow, 13 (1987) 283-294. [14] E.L. Bibeau, M. Salcudean, Subcooled void growth mechanisms and prediction at low pressure and low velocity, Int.J. Multiphase Flow, 20 (1994) 837-863. [15] O. Zeitoun, M. Shoukri, Bubble behavior and mean diameter in subcooled flow boiling, Trans. ASME, J. Heat Transfer 118 (1996) 110-116. [16] V. Prodanovic, D. Fraser, M. Salcudean, Bubble behavior in subcooled flow boiling of water at low pressures and low flow rates, Int.J. Multiphase Flow 28 (2002) 1-19. [17] I. Žun, The mechanism of bubble non-homogeneous distribution in two-phase shear flow, Nuclear Engng. Design 118 (1990) 155-162. [18] T.J. Liu, Bubble size and entrance length effects on void development in a vertical channel, Int.J. Multiphase Flow 19 (1993) 99-113. [19] A. Ohnuki, H. Akimoto, Prediction of phase distribution under bubbly flow in a large vertical pipe by multidimensional two-fluid model, Proc. 3rd Int. Conf. Multiphase Flow ICMF'98, Lyon, France, 1998. [20] T. Okawa, I. Kataoka, M. Mori, Numerical simulation of lateral phase distribution in turbulent upward bubbly two-phase flows, Nuclear Engng. Design 213 (2002) 183-197. Simulation of void fraction profile evolution in subcooled nucleate boiling 443 Strojniški vestnik - Journal of Mechanical Engineering 51(2005) 7-8, 436-444 [21] M.J. Prince, H.W. Blanch, Bubble coalescence and break-up in air-sparged bubble columns, AIChE J. 36 (1990) 1485-1499. [22] K. Sekoguchi, O. Tanaka, S. Esaki, T. Imasaka, Prediction of void fraction in subcooled and low quality boiling regions, Bull. Japan Soc. Mechanical Engineers 23 (1980) 1475-1482. [23] S.G. Kandlikar, Heat transfer characteristics in partial boiling, fully developed boiling and significant void flow regions of subcooled flow boiling, Trans. ASME, J.Heat Transfer 120 (1998) 395-401. [24] J.G.Collier, Convective Boiling and Condensation, McGraw-Hill, 1981. [25] G.Tsung-Chang, S.G. Bankoff, On the mechanism of forced-convection subcooled nucleate boiling, Trans. ASME, J.Heat Transfer 112 (1990) 213-218. [26] H.C. Unal, Maximum bubble diameter, maximum bubble-growth time and bubble-growth rate, Int.J. Heat Mass Transfer 10 (1967) 1023-1040. [27] H.J.Ivey, Relationships between bubble frequency, departure diameter and rise velocity in nucleate boiling, Int.J. Heat Mass Transfer 10 (1967) 1023-1040. [28] R. Kowe, J.C.R. Hunt, A. Hunt., B. Couet, L.J.S. Bradbury, The effects of bubbles on the volume fluxes and the pressure gradients in unsteady and non-uniform flow of liquids, Int.J. Multiphase Flow 14 (1988) 587-606. 444 Kljenak I. - Mavko B.