Nuclear Engineering and Design 237 (2007) 716–731 CFD modelling of subcooled boiling—Concept, validation and application to fuel assembly design Eckhard Krepper a,∗ , Bostjan Končar b , Yury Egorov c a Forschungszentrum Rossendorf e.V. (FZR), Institute of Safety Research, P.O.Box 510119, 01314 Dresden, Germany b Jožef Stefan Institute Jamova 39, 1000 Ljubljana, Slovenia c ANSYS, Germany GmbH, CFX-Development, Staudenfeldweg 12, D-83624 Otterﬁng, Germany Received 14 July 2006; received in revised form 19 October 2006; accepted 19 October 2006 Abstract The paper describes actual Computational Fluid Dynamics (CFD) approaches to subcooled boiling and investigates their capability to contribute to fuel assembly design. In a prototype version of the CFD code CFX a wall-boiling model is implemented based on a wall heat flux partition algorithm. It can be shown, that the wall boiling model is able to calculate the cross sectional averaged vapour volume fraction of vertical heated tubes tests with good agreement to published experimental data. The most sensitive parameters of the model are identified. Needs for more detailed experiments are established which are necessary to support further model development. The model is applied for investigation of the phenomena inside a hot channel of a fuel assembly. Here the essential phenomenon is the critical heat flux. Although subcooled boiling represents only a preliminary state towards the critical heat flux occurrence, essential parameters like swirl, cross flow between adjacent channels and concentration regions of bubbles can be determined. By calculating the temperature of the rod surface the critical regions can be identified which may later on lead to departure from nucleate boiling and possible damage of the fuel pin. The application of up-to-date CFD with a subcooled boiling model for the simulation of a hot channel enables the comparison and the evaluation of different geometrical designs of the spacer grids of a fuel rod bundle. © 2006 Elsevier B.V. All rights reserved. 1. Introduction Subcooled flow boiling occurs in many industrial applications and is characterized by large heat transfer coefficients. However, this efficient heat transfer mechanism is limited by the critical heat flux where the heat transfer coefficient decreases leading to a rapid heater temperature excursion potentially leading to heater melting and destruction. The integrity of fuel rods of a fuel assembly of a nuclear reactor core can be endangered by exceeding the critical heat flux. The critical heat flux depends on the flow parameters and can be influenced by the geometrical design of the fuel assemblies. Especially the spacer grids equipped with mixing vanes play an important role to increase the permissible heat flux. The verification of design improvements and their influence on the critical heat flux require expensive experiments. Therefore the supplementation or even the replacements ∗ Corresponding author. Tel.: +49 351 260 2067; fax: +49 351 260 2383. E-mail address: [email protected] (E. Krepper). 0029-5493/$ – see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2006.10.023 of experiments by numerical analyses are of relevant interest in fuel assembly design. In the past lots of empirical correlations for critical heat flux were developed fitted to tests and applied in purpose-developed 1D codes. But these correlations are valid only in the limited specified region of fluid conditions and fluid properties and hold for a defined geometry. Using lookup tables based on experiments, the fluid parameter validity range can be extended. But also this method is limited only to a defined geometry. The independence on the geometry can only be achieved by the application of CFD methods. The paper describes actual attempts applying CFD to solve the described task and investigates their capability to contribute to fuel assembly design. 2. Physical model Currently the most conventional CFD approach to modelling two-phase flows with significant volume fractions of both phases is the Eulerian two-fluid framework of interpenetrating continua. Phase distribution results from solving the phase- E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 Nomenclature a bubble influence factor A1F wall fraction cooled by single-phase convection A2F wall fraction cooled by quenching interfacial area density ALG CpL specific heat capacity of liquid CpG specific heat capacity of vapour dW bubble departure diameter on the wall dB mean bubble diameter in the bulk flow f bubble departure frequency g gravitational acceleration hLG interfacial heat transfer coefficient HLG evaporation enthalpy coefficient for heat transfer by quenching hQ kL liquid heat conductivity KG vapour heat conductivity active nucleation site density Na Nu Nusselt number Pr Prandtl number of liquid QE heat flux by evaporation QF heat flux by single-phase convection QQ heat flux by quenching Qtot total heat flux Re Reynolds number TL near-wall liquid temperature staturation temperature Tsat Tsub = Tsat − TL liquid subcooling Tsup = TW − Tsat superheating TW wall temperature TLW liquid characteristic near-wall temperature tW waiting time non-dimensional distance to the wall y+ UG vapour velocity UW water velocity bubble lift force Flift FWall wall lubrication force FDisp turbulent dispersion force Greek symbols α volume fraction of bubbles ρG vapour density ρL liquid density τW wall shear stress specific continuity equations for volume fractions, and a separate set of momentum equations is solved for each phase. For the steam–water bubbly flow an energy equation is solved for liquid, while vapour is assumed to be saturated everywhere. The exchange of mass, momentum and heat between phases are modelled using the correspondent source terms in the phase-specific balance equations. For the dispersed bubbly flows the interfacial momentum transfer is modelled in terms of the drag force due to the hydrodynamic resistance and the non-drag forces. The consideration of the non-drag forces namely the lift, the 717 wall lubrication and the turbulent dispersion force showed good agreement simulating air water bubbly flow (see Lucas et al., 2004). Formulation of a particular two-fluid model therefore consists of analytical or empirical correlations, used for calculating the interfacial forces, as well as the heat and mass fluxes, as the functions of the average flow parameters. Since most of these correlations are problem-specific, the range of their validity has to be kept in mind. The entire model has to be validated against experiments. A model of the internal bubbly flow with wall boiling, used in this work, is outlined in the rest of this section. 2.1. Modelling of subcooled boiling at a heated wall Subcooled boiling is observed at heated surfaces, when the heat flux applied to the wall is too high to be transferred to the core flow of liquid by the single-phase convective–conductive mechanisms. The term “subcooled” means, that the saturation temperature is exceeded only in a local vicinity of the wall, whereas the average temperature in the bulk is still below saturation. A point, where the local wall temperature reaches the saturation temperature, is considered as the onset of nucleate boiling. Steam bubbles are generated at the heated surface at nucleation sites, with the surface density of these sites depending on different factors, including the wall superheat. Further downstream the attached bubbles grow and then leave the wall at certain critical size. This critical size may depend on the surface tension and on the flow regime of the surrounding fluid. Heat transfer from the wall is then described as being carried by turbulent convection of liquid, by transient conduction due to the departing bubbles, and by evaporation. Distribution of the entire wall heat flux between these mechanisms (wall heat partitioning) can be calculated by modelling each mechanism in terms of the nucleation site density, the size of departing bubbles, their detachment frequency, and waiting time until the next bubble appears on the same site (mechanistic modelling approach). When steam bubbles move through the subcooled liquid, they condense, releasing the latent heat. The saturation temperature Tsat and the evaporation heat HLG were specified for the system pressure level and kept constant in all calculations, presented below. This assumption holds good for high pressures typical for the nuclear reactor applications. The model however can also operate with the locally calculated Tsat and HLG . For water and steam, the phase-specific densities ρL and ρG , the heat conductivity coefficients kL and kG , and the heat capacities CpL and CpG are specified at saturation conditions for the purposes of this work, while the solver supports the general temperature- and pressure-dependent properties as well. Wall and liquid temperature values participate in several constitutive model correlations in a relative form of the wall superheat above saturation Tsup =TW − Tsat , and the liquid subcooling below saturation Tsub =Tsat − TL . Here TW is the local wall temperature, and TL is the local liquid temperature of the near-wall cell. Instead of local liquid temperature TL , a characteristic temperature of liquid TLW may be used to define the magnitude of subcooling in some correlations, see later discussion. 718 E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 In technical internal flows the controlled value is the heating intensity, i.e. the wall heat flux. The given external heat flux Qtot , applied to the heated wall, is used for heating of the subcooled liquid up to and slightly above the saturation temperature (sensible heat), and for vapour generation (latent heat). The sensible heat is commonly modelled as being transported by the turbulent convection and also by the additional liquid mixing mechanism due to the bubbles, periodically emerging from the wall. The latter is termed quenching and represents the effects of transient conduction through the patches of the fresh bulk liquid, coming to the wall to replace each departing steam bubble. The quenching heat flux is noticeable on the wall area fraction A2F , influenced by steam bubbles, whereas the convection heat flux dominates on the remaining part of the wall A1F = 1 − A2F . The overall heat balance at the wall, the so-called heat partitioning, is usually written as a sum of the three parts: The following model describes the nucleate wall boiling process as the periodic release of vapour bubbles by each nucleation site: Qtot = QF + QQ + QE ṁW = ρG (1) where QF , QQ and QE are the heat flux components due to single-phase turbulent convection, quenching, and evaporation, respectively. The turbulent convection heat flux is calculated in the CFX model version using the same temperature wall function as that for the pure liquid flow without boiling (Egorov and Menter, 2004), but multiplied by the area fraction A1F : QF = A1F ρL CpL uW (TW − TL ), Ty++ (L) (2) where is analytically calculated non-dimensional temperature (Kader, 1981; Wintterle, 2004) at the non-dimensional distance from the near-wall cell y+ (L) and uW are the friction velocity. Upgrading of the boiling model with wall functions, specially customised for bubbly flows, will be the subject of further research (Končar et al., 2005b). To close the model, the remaining heat flux components QQ and QE , as well as the area fraction A2F , have to be modelled as functions of local flow parameters, including Tsup and Tsub . Also a way of estimating the local characteristic temperature of liquid TLW in Tsub must be defined. After that the Eq. (1) can be solved iteratively for the local wall temperature TW , which satisfies the wall heat flux balance. Similar to the model by Kurul and Podowski (1991), the CFX model represents QQ in terms of the quenching heat transfer coefficient hQ (3) To obtain the grid independent solution for QQ and taking into account the self-similarity of non-dimensional temperature profiles at different y+ , the characteristic temperature of liquid TLW in Eq. (3) can be calculated as: TLW = TW − Ty++ (LW) Ty++ (L) (TW − TL ). (4) The evaporation heat flux QE is obtained via the evaporation mass flux on the wall QE = ṁW HLG . (6) A2F (7) Here dW is the bubble departure diameter, Na the surface density of active nucleation sites, f the departure frequency, and a is the bubble influence factor. A value of 2 is taken for the bubble influence factor a, participating in the area fraction A2F (7). The effect of area overlapping between the neighbouring sites is here not modelled, therefore the evaporation rate ṁW is calculated using the area-limited variant of the relation (6): 2 dW A2F Na f 3 a2 (5) (8) In this work basically the same correlations for dW , Na , f, a, and hQ are applied as those used by Kurul and Podowski (1991). The active nucleation site density Na depends on the wall surface roughness, on the wettability of the solid-liquid pair and above all on the wall superheat. For water and steel Na is correlated according to Lemmert and Chawla (1977) depending solely on the wall superheat as Na [m−2 ] = (185 Tsup [K])1.805 Ty++ (L) QQ = A2F hQ (TW − TLW ) 3 πdW Na f 6 (adW )2 = min π ,1 4 ṁW = ρG (9) For the bubble departure diameter dW , a following correlation was applied (Tolubinsky and Kostanchuk, 1970): −Tsub,LW , 1.4 [mm] (10) dW = min 0.6 [mm] exp 45 K The model is derived for the high-pressure water boiling experimental data with the upper limit for the bubble departure diameter (dBW = 1.4 mm). The correlation (10) depends solely on characteristic liquid subcooling, which makes this correlation case dependent. In the literature, also mechanistic models based on the microscopic phenomena may be found. The model of (Ünal, 1976) presumes that the bubble is generated by evaporation of the superheated liquid micro-layer, while at the same time the top of the bubble is subjected to condensation. The bubble may slide or collapse on the heated surface during its lifetime. After reaching the maximum size, the bubble departs from the wall. Further investigation is necessary, since the model of bubble departure diameter dW has a strong influence on the share of the heat transfer components and also on the calculated amount of steam production. A simple estimation of the bubble departure frequency as the terminal rise velocity over the bubble departure diameter is adopted here (Cole, 1960): 4g(ρL − ρG ) (11) f = 3dW ρL The quenching heat transfer coefficient is calculated using the analytical solution for one-dimensional transient conduction, as E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 suggested by Mikic and Rohsenow (1969): 2 hQ = √ f tW kL ρL CpL π (12) where tW is waiting time between the bubble departure and the appearance of a next bubble at the same nucleation site. A simple assumption by Kurul and Podowski (1991) is adopted also here, where the waiting time takes 80% of the bubble departure period: tW = 0.8 f (13) Definition of the characteristic temperature TLW is a key point in the CFD-implementation of the wall boiling model. Its importance arises from the spatially averaged representation of the boiling phenomena by the mechanistic model, which must be combined with the detailed local nature of the CFD method. The constitutive correlations for the quenching heat flux (Eqs. (3) and (12)) and for the bubble departure diameter (Eq. (10)) were originally formulated for the one-dimensional thermo-hydraulic models in terms of the mean-flow liquid temperature. Some of the other popular correlations, not used in this work, also depend on the core flow velocity. Formal implementation of such correlations as the CFD wall boundary conditions assumes replacing the non-local core flow quantities by the near-wall local ones. It can only be done for the extremely coarse grids, with the first near-wall grid cell covering the whole boundary layer thick- 719 ness. With the near-wall grid, which is adequate for the CFD turbulent wall treatment, such a straightforward approach can significantly overestimate the vapour generation, since the CFD solution attempts to resolve the internal superheated part of a temperature boundary layer in the liquid phase (Menter, 2003). In CFX, the characteristic temperature TLW is taken from the reconstructed temperature profile using the temperature wall function (4) and the given value of the non-dimensional distance y+ to the wall. In this work a constant value of 250 is used for y+ , although more sophisticated choices are also possible, for example calculating y+ on the base of the local bubble departure diameter, or on the base of the superheated liquid layer thickness. The model, implemented in CFX, is found virtually invariant to the near-wall grid cell size in a wide range covering the logarithmic layer. Besides the improved heat and mass transfer representation, this implementation enables the use of the adequate grids for the velocity wall function, which is a prerequisite for the accurate calculation of the wall shear stress and the pressure drop in internal flows. Steam is assumed to be at saturation condition. Within the subcooled liquid (TL < Tsat ) steam is condensing with the mass transfer rate per unit volume: ṁ = max hLG (Tsat − TL )ALG ,0 HLG (14) Fig. 1. Calculated distributions of void fraction (a) and temperature (b) in the channel; (c) measured vs. calculated averaged void fraction, averaged temperature, axis temperature and wall temperature. 720 E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 With superheated liquid, fluid is evaporating at the rate: ṁ = max hLG (TL − Tsat )ALG ,0 HLG 2.2. Modelling of the momentum transfer (15) ALG is the interfacial area, and hLG is the interfacial heat transfer coefficient, calculated according to Ranz and Marshall (1952): hLG = kL kL Nu = (2 + 0.6Re1/2 Pr1/3 ) dB dB (16) To close the phase transition model in the bulk bubbly flow with a mean bubble diameter dB , Kurul and Podowski (1991) and also Anglart et al. (1997) proposed to calculate the bubble diameter dB locally as a linear function of liquid subcooling Tsub : dB = dB1 (Tsub − Tsub,2 ) + dB2 (Tsub,1 − Tsub ) Tsub,1 − Tsub,2 (17) For typical nuclear energy applications these authors proposed reference bubble diameters at the two reference subcooling conditions: dB1 = 0.1 mm at Tsub,1 = 13.5 K and dB2 = 2 mm at Tsub,2 = 5 K. The bubble size in the bulk has a direct influence on the interfacial area density and on the condensation respective evaporation rate in the bulk. For the liquid phase the Shear Stress Transport (SST) turbulence model by Menter (1994) was applied, which operates for the considered here internal flows in long pipes and channels similar to the k–ε model. The influence of the gas bubbles on turbulence in liquid was modelled using Sato et al.’s (1981) eddy viscosity model for bubble-induced turbulence. Besides the drag forces, representing the flow resistance and modelled here using the correlation by Ishii and Zuber (1979), the non-drag forces also need to be modelled to predict the non-homogeneous flow structure. Namely the turbulent dispersion, the lift force, and the wall lubrication force should be considered. In the following expressions the forces acting on gas bubbles are given. The non-drag forces were investigated analysing lots of experiments of adiabatic air/water bubbly flow in a vertical tube for a large parameter range of superficial velocities performed in FZ-Rossendorf (see Lucas et al., 2005). The applied turbulent dispersion force is based on the Favre average of the interfacial drag force (Burns et al., 2004) and is calculated for a two phase flow according to FDisp = − 3CD μt grad α (UG − UL ) 4dB σt 1−α Fig. 2. Calculated heat flux partition (a) and bubble size diameter at departure dW and in the bulk dB (b). (18) E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 721 Fig. 3. Condensation rate (a) and cross sectional averaged condensation, bulk evaporation rate and interfacial area (b). with CD as the bubble drag coefficient, μt the dynamic eddy viscosity of liquid, and σ t the turbulent Schmidt number for the volume fraction of the dispersed phase. The lift force is defined in a standard way as: using the gas and liquid velocities UG and UL and the liquid density ρL . Tomiyama et al. (1995), Tomiyama (1998) has performed extensive investigations of the lift force coefficient Clift as a function of the bubble size in the air–water two-phase flow. He has found a sign change of Clift : when growing dB passes a critical value, the lift coefficient becomes negative. For air water flow at atmospheric pressure this value amounts to 5.8 mm. The evaluation of Tomiyamas correlation for Clift (Tomiyama, 1998) and experiments performed at TOPFLOW (Prasser et al., 2005) show that for steam water flow at higher pressure the critical bubble size is shifted to smaller values: At 5 MPa to about 3.5 mm and at 15 MPa to about 2 mm. Nevertheless, the bubble size in Fig. 4. Measured (symbols) and calculated (solid lines) vapour volume fractions dependent on the quality x. Fig. 5. D = 12 mm, Q = 1.0 MW/m2 , G = 1000 kg/(m2 s). FLift = −Clift ρL (UG − UL ) × rot UL (19) 722 E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 and another correlation by Tomiyama (1998) FWall T Fig. 6. D = 12 mm, Q = 0.8 MW/m2 , P = 7 MPa. the here presented application cases, which amounts at departure in the order of 0.5 mm and in the bulk in the order of 0.2–1.5 mm will be lower than the critical value. Investigations with air water flow indicate, that the system of non-drag forces, which yields very good results with velocities of about 1 m/s overestimates the gas fraction wall peak at high water velocities >4 m/s (see Lucas et al., 2004). The wall lubrication force pushes the bubbles away from the wall. The available model correlations have been tested only for adiabatic air water bubbly flows without bubble generation at the wall, and therefore their use with the wall boiling model is debatable. Known correlations are a correlation by Antal et al. (1991) FWall A = − ρW α dB CW1 − CW2 dB y 2 n Urel Fig. 7. D = 12 mm, P = 7 MPa, G = 1000 kg/(m2 s). (20) dB α = −CW 2 1 1 − y2 (D − y)2 2 ρL Urel n (21) with n as the normal vector to the wall, dB as the bubble size diameter, Urel as the velocity difference between the phases, and α as the gas volume fraction. The force is dependent on the distance to the wall y. Antal et al. (1991) proposed an inverse proportionality to y, whereas Tomiyama (1998) proposed an inverse quadratic proportionality to y. The latter model has been extensively validated for adiabatic air water bubbly flows, whereas other investigations showed, that the model of Antal seems to be more appropriate for wall boiling flows (see Končar et al., 2005a). A final statement on this issue cannot yet be done and further investigations are necessary. Although the implementation of the wall lubrication force is necessary for the adiabatic two-phase flows, as it reproduces the void fraction peak near the wall (Lucas et al., 2004) its use at high-pressure wall boiling conditions may be questionable. Pierre and Bankoff (1966) have found in their experiments no void peak near the wall. Also, no void peak has been observed at boiling of refrigerants, where the fluid properties are similar to boiling of water at high pressures (Roy et al., 2002). At wall boiling the spherical-cap bubbles continuously grow and/or slide on the heated wall. According to its formulation a zero void fraction at the wall is required. The Antal or Tomiyama formulation of the wall force (Eqs. (20) and (21)) formally contradict with the non-zero wall source of dispersed phase, since their integrals across the boundary layer are infinitely large. Still, the use of wall lubrication force (Končar et al., 2004) was proved to be advantageous for the conditions of low-pressure boiling of water, where the void fraction peak was observed (Bartel et al., 2001). Further research is necessary to improve the existing wall force models. Fig. 8. D = 12 mm, P = 15 MPa, G = 2000 kg/(m2 s). E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 3. Validation of the model An experimental arrangement suited for the comparison of measured and calculated results is an upward water flow through a heated tube. The boundary conditions can be clearly specified and the resulting cross sectional averaged steam volume fraction can be measured by gamma densitometry with acceptable accuracy. Numerous experiments of this kind are reported in the literature. Examples are the experiments published by Bartolomej and Chanturiya (1967) and Bartolomej et al. (1982). The boiling model, described in the previous chapter, was applied to simulate these experiments to show the plausibility of the calculated results. To identify the range of validity of applied boiling model the calculated cross sectional averaged steam volume fractions were compared with the large set of measurements at different operating conditions 3.1. Benchmark case For the benchmark calculation the Bartolomej experiment performed in 2 m long heated tube with the inner diameter of 15.4 mm was selected. The heat flux was 5.7 × 105 W/m2 and the mass flow of the water at the pressure of 4.5 MPa amounted to 900.0 kg/(s m2 ). In the calculations the inlet subcooling was set to 58.2 K so that the thermodynamic quality x = 0 723 was reached 1.75 m downstream of the tube inlet. For this particular experiment also the wall temperature along the tube was published. The tube was modelled in 2D cylindrical geometry with 20 grid cells in radial and 150 grid cells in axial direction. The grid refinement analysis performed in our previous work (Končar et al., 2005a) showed that this grid is a good compromise between the numerical accuracy and the computational effort. The following closure models were used for the reference calculation: the SST turbulence model for the liquid phase, the Sato model for the bubble induced turbulence, the Ishii–Zuber model of the interfacial drag, the Tolubinsky model for the bubble departure diameter (Eq. (10)) and the model for bulk bubble diameter size according to Eq. (17). The Favre-averaged turbulent dispersion force by Burns et al. (2004) was the only considered non-drag force. Since for high-pressure boiling of water there is no available measurements of local two-phase flow parameters, only the simplest closure models have been selected here. The decision on the selected set of correlations is also based on good agreement with experiments over a wide range of flow conditions, as shown in the next section. Fig. 1a shows the calculated distribution of the steam volume fraction in the bulk and Fig. 1b of the water temperature. The presentation is stretched in radial direction. The tube is heated from the right side and the central symmetry axis is on the left Fig. 9. Influence of doubling the bubble diameter at departure on the wall heat flux partition (a) and on the condensation rate (b) (dotted lines—reference case). 724 E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 Fig. 10. Influence of doubling the bubble diameter in the bulk on the condensation rate (a) and on the vapour volume fraction (b) (dotted lines—reference case). side. Fig. 1c shows the comparison of measured (symbols) and calculated (solid lines) cross sectional averaged values of steam volume fraction and water temperature and also the wall temperature values along the axial tube length. Both, the point of boiling inception and the increasing of the vapour volume fraction over the height show good agreement to the measured data. Also the calculated temperatures correspond very well to the experiment. To show the plausibility of the modelled results, further parameters are extracted from the solution. The calculated wall heat partitioning (Fig. 2a) shows, that at the beginning almost all heat flux is transferred by the single-phase convection, whereas Fig. 11. (a) Radial distribution considering all non-drag forces. The dotted line shows the reference case, where wall and lift bubble force were not considered. (b) Influence of the non-drag forces on the cross sectional averaged vapour volume fraction (dotted lines—reference case). E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 Fig. 12. Vapour volume fraction distribution in the cross section at z = 0.475 m without consideration of an inlet swirl. with increasing tube length quenching and evaporation become larger. The cross sectional averaged values of bubble departure diameter and local bubble diameter in the bulk are compared in Fig. 2b. According to Eq. (17) the bulk bubble diameter strongly depends on the local water temperature in the bulk therefore it increases significantly towards the end of the channel where bulk evaporation occurs. Fig. 3a shows the condensation rate calculated according to Eq. (14), which depends on the liquid subcooling, the vapour volume fraction and the interfacial area. Fig. 3b shows cross sectional averaged values of condensation rate, bulk evaporation rate and interfacial area along the channel height. 725 presented model setup. The calculations were performed with the same set of models as for the benchmark case. Fig. 4 shows measured (symbols) and calculated (solid lines) cross sectional averaged steam volume fractions dependent on the thermodynamic quality x for different pressures (see Bartolomej and Chanturiya, 1967). The first case at 1.5 MPa corresponds to the test presented in the previous chapter. The tube diameter was 15.4 mm, the heat flux 3.8 × 105 W/m2 and the mass flow amounted to 900.0 kg/(s m2 ). The point of boiling inception is very well predicted and also the trends of the volume fraction development with increased thermodynamic quality show good agreement with experiments. Figs. 5–8 show the results of further test series performed in a tube with inner diameter of 12 mm (Bartolomej et al., 1982). The pressure effect is presented in Fig. 5. A reasonable agreement with experiments was obtained at pressures between 3 and 11 MPa, whereas the averaged void fraction is somewhat underpredicted at 15 MPa. The influence of mass flux at constant pressure and heat flux is presented in Fig. 6. A significant under-estimation of void fraction was found at the lowest mass flux of G = 400 kg/(m2 s), whereas a good agreement with experiments may be observed at higher mass fluxes. The influence of heat flux at constant values pressure p = 7 MPa and mass flux G = 1000 kg/(m2 s) is presented in Fig. 7. A very good prediction may be observed at moderate heat fluxes, whereas some discrepancy occurs at higher heat flux values of 1.7 and 2 MW/m2 . 3.2. Comparison with experiments—range of validity To investigate the range of validity of the subcooled boiling model at different pressure, mass flow and heat flux conditions, numerous experiments performed in different circular tubes and published by Bartolomej (1967, 1982) were compared to the Fig. 13. Velocity field at the inlet: axial velocity 5 m/s, velocity at the circumference of the swirl 3 m/s (case d). Fig. 14. Axial profiles of the swirl (see Eq. (21)) for the investigated cases (vz = 5m/s). 726 E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 Fig. 15. Stream line presentation of the simulated heat channel. The heat flux effect at higher pressure p = 15 MPa and mass flux G = 2000 kg/(m2 s) may be observed in Fig. 8. Here the void fraction is over-predicted at most of the heat flux values, whereas the discrepancy is the smallest at the highest heat flux (2.2 MW/m2 ). Discrepancies at high pressures (15 MPa) may be caused due the use of simple model for bubble departure diameter (Tolubinsky and Kostanchuk, 1970), which is not pressure dependent. Nevertheless, in all presented cases, the qualitative trends are reasonably predicted. Summarizing this comparison, the model presented in the previous chapter is best suited for pressures from 3 to 11 MPa, for heat fluxes up to 1.2 MW/m2 and for mass flow rates at about 1000 kg/(m2 s). the cross sectional averaged void fraction however is partially compensated by a higher condensation rate (see Fig. 9b). The bubble size in the bulk has via the interfacial area an influence on the condensation rate (see Fig. 10a). The effect on the vapour volume fraction (see Fig. 10b) is partially compensated by the radial temperature profile (see Fig. 11). Most of the vapour is found in areas of temperatures near the saturation temperature and small condensation rates. Although the influence of the bubble sizes on the vapour volume fraction is limited, the calculations show the influence 3.3. Identiﬁcation of most sensitive model parameters In the actual calculations a quite simple correlation for the bubble size diameter at departure dW fitted to experiments was used (see Eq. (10)). In the literature correlations based on mechanistic models of the microscopic phenomena are known (e.g. Ünal, 1976) and their application may extend the parameter range of validity. To demonstrate the sensitivity of the calculated vapour volume fraction on the modeling of dW , the bubble bubble departure diameter was doubled. The bubbles grow longer at the wall before they leave into the bulk. The calculation with the doubled bubble size in the Figs. 9 and 10 is represented by the solid lines, whereas the dotted lines show the reference case. In the case of doubled diameter a higher fraction of the heat flux contributes to evaporation (see Fig. 9a). The effect on Fig. 16. Radial profiles of the non-drag forces for the rod bundle geometry at VR = 3m/s. E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 727 on the relation of the microscopic mechanisms. The bubble size at departure has a large influence on the wall heat partition. The simulated bubble size in the bulk influences the condensation rate mainly via the interfacial area density. Unfortunately the measurement of the bubble sizes by actual techniques is very difficult. Difficulties arise due to very narrow channels of about 10 mm, where the size of the vapour bubbles is expected in the range between 0.5 and 2 mm. Furthermore, the pressure and temperature of the working fluid typical for nuclear reactor is very high. The design of experiments using model fluids enables the enlargement of the geometry and easier determination of cross sectional vapour distribution. Recent attempts applying tomographic methods are very promising. 3.4. Inﬂuence of non-drag forces To show the influence of the lift and the wall force, corresponding simulations were performed. First test calculations revealed, that the assumed lift coefficient Clift = 0.288 like in air water adiabatic flow is too high. Caused by the shifted water velocity profile by the generated steam all generated vapour remained at the rod surface. Therefore, for the calculations presented in this study, a smaller lift force coefficient of Clift = 0.06 was assumed. A wall force according to Eq. (20) (Antal et al., 1991) with the coefficients CW1 = −0.025 and CW2 = 0.075 was considered. The non-drag forces have influence on the radial vapour volume fraction profiles. In Fig. 11a radial profiles of the liquid velocity (blue line, left axis), gas velocity (green line, left axis), gas volume fraction (red line, right axis) and liquid temperature (black line, left temperature axis) are presented. At z = 0.5 m the flow is still single phase and a typical profile of a liquid flow in a tube can be found. At 1.75 m steam is present. The vapour velocity flow profile is shown by the green line, the liquid profile is seen accelerated by the vapour. The consideration of the bubble forces in the calculation caused that the maximum value of the vapour volume fraction decreases and is shifted somewhat away from the wall. The maximum liquid temperature (black solid line) is calculated just at the wall. The influence on the cross sectional averaged vapour volume fraction is almost negligible (see Fig. 11b). Fig. 17. Vapour volume fraction distribution in the cross section at z = 0.475 m considering a swirl (case d). of 12.6 mm. The thermal hydraulic and transport water properties were set for a pressure of 15.7 MPa, typical for PWR conditions. The heat flux at the rod surface was assumed to be 1.0 × 106 W/m2 and the subcooling at the inlet was set to 12 K expecting the generation of vapour in the simulated section. The axial water velocity was set to Uz = 5 m/s. The faces at the low and high x respective at low and high y were simulated as periodic boundary conditions. 4. Simulation of a hot channel of a fuel assembly In this chapter CFX calculations of a hot channel between four rods are presented. Although the CFD modelling of critical heat flux is not yet possible, the simulations shall demonstrate the capability of CFD supporting the fuel assembly design. In the calculations only subcooled boiling is simulated, which is here considered as a preliminary phenomenon towards departure of nucleate boiling (DNB). DNB might occur at the thermal hydraulic conditions of a PWR. A situation was investigated, when at full power and full pressure the inlet temperature rise caused undesired boiling in the channel. A section between two spacer grids having a length of z = 0.5 m was simulated. The grid represents a subchannel between 4 rods having a diameter of 9 mm and a rod distance Fig. 18. Influence of the swirl on the cross sectional averaged vapour volume fraction and water temperature for case d). 728 E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 The same model setup for mass, momentum and heat transfer as in the benchmark case of Section 3.1 was applied. The result is a vapour volume fraction distribution in the channel cross section at z = 0.475 m (see Fig. 12). The figure shows maximum vapour generation at the smallest distance between the rods. To investigate the effect of a swirl caused by mixing vanes, additional horizontal velocity components UX and UY were given at the inlet according to Eq. (22). y π = −UR Ux = UR cos α + 2 r x π Uy = UR sin α + = UR , 2 r for 1 mm < r < 3.8 mm (22) Besides the reference case (a) with zero transversal velocity (UR = 0 m/s), the following cases with different transversal velocities at the circumference of the circular swirl VR were simulated: (b) 1 m/s, (c) 2 m/s and (d) 3 m/s to simulate examples for different geometrical design of mixing vanes. Fig. 13 shows the resulting inlet velocity field for the case (d) with UR = 3 m/s. Fig. 14 presents the relation of the horizontal velocity components to the overall value of the velocity corresponding to Eq. (23) along the vertical axis z. Ux2 + Uy2 swirl = Ux2 + Uy2 + Uz2 (23) Due to the centrifugal force, the heavier fluid component – the water – is pushed outwards, whereas a large amount of the lighter component – the vapour – is accumulated in the centre of the channel. The streamline presentations of Fig. 15a for water and Fig. 15b for vapour clearly show this phenomenon. The centrifugal force Fcent attracting at the gaseous phase can be calculated according to ×ω F cent = −ρL U (24) =∇ ×U (25) and the vorticity ω. with the velocity vector U Fig. 16 shows the radial profiles of the three bubble forces √ projected to a normal vector (1/ 2) (1; 1; 0) pointing from the centre of the channel to the centre of a rod. FDisp is calculated according to Eq. (18) and Flift according to Eq. (19). The figure shows that in this geometry and under these flow conditions the flow is mainly determined by the centrifugal force and the turbulent dispersion force. The latter force acts mainly in the Fig. 19. Cross sectional averaged heat flux components (a) and superheating temperature at the rod surface (b) for the case d) (solid lines) and the reference case without swirl (dotted lines). E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 centre of the channel. Here a large vapour volume fraction and a large gradient of the vapour volume fraction can be found. The lift force plays only a limited role. Fig. 17 shows the calculated vapour volume fraction distribution at z = 0.475 m considering the counter clockwise swirl. In the shadow of the rods bubble accumulations are calculated. In the calculations like in the basic case only a monodispersed bubble size distribution in the bulk dependent on the liquid temperature according Eq. (17) was modelled. Further work has to 729 consider however for this high vapour volume fraction in the centre bubble coalescence and the occurrence of larger bubbles. Hence in the actual calculation the interfacial area and consequently the condensation in the bulk are overestimated. In the reality, higher vapour volume fractions might be expected. In Fig. 18 the influence of the swirl on the axial profile of the cross sectional averaged vapour volume fraction for the case (d) is presented. Indeed the cross sectional averaged vapour volume fraction is reduced by the swirl as desired—in this example Fig. 20. Distribution of the wall superheating temperature at the surface of rod 1 without swirl (a) and considering a swirl (cases b–d) (for the definition of the angle Φ, here in radians, see Fig. 11). 730 E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 from about 8% to about 6.7%, whereas the averaged water temperature in the channel is almost not influenced. The tangential velocity near the wall is increased by the swirl. As a consequence, higher portions of the wall heat flux are transformed by single-phase convection whereas the share of evaporation is decreased (see Fig. 19a). Furthermore the cross sectional averaged wall superheating temperature is reduced (see Fig. 19b). In the presented example the reduction at z = 0.1 m amounts to about 1.5 K. Very informative is the distribution of the superheating Tsup at the rod surface of rod 1 shown in Fig. 20 for the four cases (a–d). Without swirl, a symmetric distribution is found showing maximum values at locations with the smallest distance of the adjacent rods. In contrast, the assumed swirl generates a non-equally distributed wall superheating. Hot spots (e.g. at z = 0.15, Φ = 0.3) might be the first locations, where later on the critical heat flux could occur. 5. Summary and conclusions Capabilities of the up-to-date two-fluid CFD method of simulating steam–water bubbly flows with phase transition are demonstrated. To simulate boiling at a heated wall, the phenomena on the micro-scale have been modelled by appropriate closure relations. The given heat is removed from the wall by different mechanisms. The size of bubbles, emerging from the wall, has been proven to have a large influence on the heat distribution between the different heat flux components. The size of bubbles in the bulk is correlated to the local liquid subcooling, since it determines the interfacial area and the condensation rate of the vapour bubbles. Further improvements of the microscopic models are necessary. The model for the bubble size at departure has to be replaced by a mechanistic model. The simulation of the bubble size in the bulk has to consider bubble coalescence and fragmentation. The experimental design needs to be optimised to deliver the relevant information for the CFD model improvement. The determination of the spatial vapour volume distribution is necessary to support further model developments, especially modelling of the momentum exchange between gas bubbles and the liquid. The determination of the bubble size would be desirable for correct simulation of the bubble behaviour at the heated surface and in the bulk. Nevertheless, despite its current developmental state, the used CFD model yields to correct quantitative results and predicts right trends in parametric studies. Although only a preliminary state towards the critical heat flux is simulated, essential parameters like the swirl, the cross flow between adjacent channels and concentration regions of bubbles can be determined. By predicting the temperature at the rod surface, critical regions can be identified, which might later on lead to departure from nucleate boiling and possible damage of the fuel pin. The avoiding of uneven temperature distributions at the rod surface can be considered as a design criterion for the spacer grids. The application of actual CFD with a subcooled boiling model for the simulation of a hot channel enables at least the comparison and the evaluation of different geometrical designs of the spacer grids. References Anglart, H., Nylund, O., Kurul, N., Podowski, M.Z., 1997. CFD prediction of flow and phase distribution in fuel assemblies with spacers. NURETH-7, Saratoga Springs, New York, 1995. Nucl. Eng. Des. 177, 215–228. Antal, S.P., Lahey, R.T., Flaherty, J.E., 1991. Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Int. J. Multiphase Flow 7, 635–652. Bartel, M.D., Ishii, M., Masukawa, T., Mi, Y., Situ, R., 2001. Interfacial area measurements in subccoled flow boiling. Nucl. Eng. Des. 210, 135–155. Bartolomej, G.G., Chanturiya, V.M., 1967. Experimental study of true void fraction when boiling subcooled water in vertical tubes. Thermal Engineering, vol. 14, pp. 123–128 (translated from Teploenergetika, no. 2, vol. 14, 1967, pp. 80–83). Bartolomej, G.G., et al., 1982. An experimental investigation of true volumetric vapour content with subcooled boiling in tubes. Thermal Engineering, vol. 29, pp. 132–135 (translated from Teploenergetika, no. 3, vol. 29, 1982, pp. 20–23). Burns, A.D., Frank, T., Hamill, I., Shi, J.-M., 2004. The favre averaged drag model for turbulence dispersion in Eulerian multi-phase flows. In: 5th Int. Conf. on Multiphase Flow, ICMF’2004. Yokohama, Japan. Cole, R., 1960. A photographic study of pool boiling in the region of the critical heat flux. AIChE J. 6, 533–542. Egorov, Y., Menter, F., 2004. Experimental Implementation of the RPI Wall Boiling Model in CFX-5.6. Technical Report ANSYS/TR-04-10. ANSYS Germany GmbH. Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25, 843–855. Kader, B.A., 1981. Temperature and concentration profiles in fully turbulent boundary layers. Int. J. Heat Mass Transfer 24, 1541–1544. Končar, B., Kljenak, I., Mavko, B., 2004. Modelling of local two-phase flow parameters in upward subcooled flow boiling at low pressure. Int. J. Heat Mass Transfer 47, 1499–1513. Končar, B., Krepper, E., Egorov, Y., 2005a. CFD Modelling of Subcooled Flow Boiling for Nuclear Engineering Applications International Conference Nuclear Energy for New Europe, Bled, Slovenia, September 5–8, 2005 (Paper 140). Končar, B., Mavko, B., Hassan, Y.A., 2005b. Two-phase wall function for modeling of turbulent boundary layer in subcooled boiling flow. In: 11th International Topical Meeting on Nuclear Reactor Thermal-Hydraulics (NURETH-11), Avignon, France, October 2005 (Paper 443). Kurul, N., Podowski, M.Z., 1991. On the modeling of multidimensional effects in boiling channels. In: ANS Proceedings of 27th National Heat Transfer Conference, Minneapolis, MN. Lemmert, M., Chawla, J.M., 1977. Influence of flow velocity on surface boiling heat transfer coefficient. In: Hahne, E., Grigull, U. (Eds.), Heat Transfer in Boiling. Academic Press and Hemisphere, ISBN 0-12-314450-7, pp. 237–247. Lucas, D., Shi, J.-M., Krepper, E., Prasser, H.-M., 2004. Models for the forces acting on bubbles in comparison with experimental data for vertical pipe flow. In: 3rd International Symposium on Two-Phase Flow Modelling and Experimentation, Pisa, Italy (Paper ha04). Lucas, D., Krepper, E., Prasser, H.-M., 2005. Development of co-current airwater flow in a vertical pipe. Int. J. Multiphase Flow 31, 1304–1328. Menter, F.R., 2003. Personal communication. Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32, 8. Mikic, B.B., Rohsenow, W.M., 1969. A new correlation of pool-boiling data including the fact of heating surface characteristics. ASME J. Heat Trans. 91, 245–250. Pierre, C.C., Bankoff, S.G., 1966. Vapor volume profiles in developing twophase flow. Int. J. Heat Mass Transfer 10, 237–249. Prasser, H.-M., Beyer, M., Carl, H., Gregor, S., Lucas, D., Pietruske, H., Schütz, P., Weiss, F.-P., (Paper 399) 2005. Evolution of the structure of a gas- E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731 liquid two-phase flow in a large vertical pipe. In: NURETH-11, Avignon, France. Ranz, W.E., Marshall, W.R., 1952. Evaporation from drops. Chem. Eng. Progr. 48 (3), 141–146. Roy, R.P., Kang, S., Zarate, J.A., Laporta, A., 2002. Turbulent subcooled boiling flow—experiments and simulations. J. Heat Trans. 124, 73–93. Sato, Y., Sadatomi, M., Sekoguchi, K., 1981. Momentum and heat transfer in two-phase bubble flow. Int. J. Multiphase Flow 7, 167–177. Tolubinsky, V.I., Kostanchuk, D.M., 1970. Vapour bubbles groth rate and heat transfer intensity at subcooled water boiling; Heat Transfer 1970, Preprints of papers presented at the 4th International Heat Transfer Conference, vol. 5, Paris (Paper No. B-2.8). 731 Tomiyama, A., et al., 1995. Effects of Eötvös number and dimensionless liquid volumetric flux on lateral motion of a bubble in a laminar duct flow. Adv. Multiphase Flow, 3–15. Tomiyama, A., 1998. Struggle with computational bubble dynamics. In: Third International Conference on Multiphase Flow, ICMF’98, Lyon, France, June 8–12, 1998. Ünal, H.C., 1976. Maximum bubble diameter, maximum bubble growth time and bubble growth rate. Int. J. Heat Mass Trans. 19, 643–649. Wintterle, T. 2004. Development of a numerical boundary condition for the simulation of nucleate boiling at heated walls. Diploma Thesis University. Stuttgart (IKE-8-D-014).