CFD modelling of polydispersed bubbly two-phase flow around an obstacle

Nuclear Engineering and Design 239 (2009) 2372–2381

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

CFD modelling of polydispersed bubbly two-phase flow around an obstacle Eckhard Krepper a,∗ , Matthias Beyer a , Thomas Frank b , Dirk Lucas a , Horst-Michael Prasser c a

Forschungszentrum Dresden-Rossendorf e.V., (FZD), Institute of Safety Research, P.O. Box 510119, 01314 Dresden, Germany ANSYS Germany GmbH Staudenfeldweg 12, 83624 Otterfing, Germany c ETH Zürich, Institut für Energietechnik, Sonneggstrasse, 8092 Zürich, Switzerland b

a r t i c l e

i n f o

Article history: Received 20 July 2008 Received in revised form 26 May 2009 Accepted 9 June 2009

a b s t r a c t A population balance model (the Inhomogeneous MUSIG model) has recently been developed in close cooperation between ANSYS-CFX and Forschungszentrum Dresden-Rossendorf and implemented into the CFD-Code CFX [Krepper, E., Lucas, D., Prasser, H.-M, 2005. On the modelling of bubbly flow in vertical pipes. Nucl. Eng. Des. 235, 597–611; Frank, T., Zwart, P.J., Shi, J.-M., Krepper, E., Rohde, U., 2005. Inhomogeneous MUSIG Model—a population balance approach for polydispersed bubbly flows, International Conference “Nuclear Energy for New Europe 2005”, Bled, Slovenia, September 5–8, 2005; Krepper, E., Beyer, M., Frank, Th., Lucas, D., Prasser, H.-M., 2007. Application of a population balance approach for polydispersed bubbly flows, 6th Int. Conf. on Multiphase Flow Leipzig 2007, (paper 378)]. The current paper presents a brief description of the model principles. The capabilities of this model are discussed via the example of a bubbly flow around a half-moon shaped obstacle arranged in a 200 mm pipe. In applying the Inhomogeneous MUSIG approach, a deeper understanding of the flow structures is possible and the model allows effects of polydispersion to be investigated. For the complex flow around the obstacle, the general structure of the flow was well reproduced in the simulations. This test case demonstrates the complicated interplay between size dependent bubble migration and the effects of bubble coalescence and breakup on real flows. The closure models that characterize the bubble forces responsible for the simulation of bubble migration show agreement with the experimental observations. However, clear deviations occur for bubble coalescence and fragmentation. The models applied here, which describe bubble fragmentation and coalescence could be proved as a weakness in the validity of numerous CFD analyses of vertical upward two-phase pipe flow. Further work on this topic is under way. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Many flow regimes in Nuclear Reactor Safety Research are characterized by multiphase flows, with one phase being a liquid and the other phase consisting of gas or vapour of the liquid phase. The flow regimes found in vertical pipes are dependent on the void fraction of the gaseous phase, which varies with increasing void fraction from bubbly flow to slug flow, churn turbulent flow, annular flow and finally to droplet flow at highest void fractions. In the regimes of bubbly and slug flows, a spectrum of different bubble sizes is observed. While dispersed bubbly flows with low gas volume fraction are mostly mono-dispersed, an increase in the gas volume fraction leads to a broader bubble size distribution due to breakup and coalescence of bubbles. The interfacial area is crucially important for the transfer of mass, momentum or heat between the continuous and the dispersed phase; thus, the interfacial area

∗ Corresponding author. Tel.: +9 351 260 2067; fax: +49 351 260 12067. E-mail addresses: [email protected] (E. Krepper), [email protected] (T. Frank), [email protected] (D. Lucas), [email protected] (H.-M. Prasser). 0029-5493/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.nucengdes.2009.06.015

cannot be simply modelled under the assumptions of just the void fraction and a mean diameter. Moreover, the forces acting on the bubbles may depend on their individual size. This is the case not just for drag, but also for non-drag forces. Among the forces leading to lateral migration of the bubbles, i.e. acting in perpendicular direction with respect to the main drag force, the bubble lift force was found to change the sign as the bubble size varies. Consequently, in the context of pipe flows, this leads to a radial separation between small and large bubbles and to enhanced coalescence of large bubbles, which migrate towards the pipe center and form even larger Taylor bubbles or slugs. An adequate modelling approach must consider all these phenomena. A generalized Inhomogeneous Multiple Size Group (MUSIG) model was applied here, where the dispersed gaseous phase is divided into N inhomogeneous velocity groups (phases), each of these groups is then subdivided into Mj bubble size classes (see Frank et al., 2006, 2007). The bubble breakup and coalescence processes are taken into account by appropriate models. The model was developed, tested and first validated via comparisons with experimental data obtained for vertical pipe flow in the cases of air–water and steam–water at saturation conditions

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

Nomenclature ai CL d Eo f F˙ c,b FB/C FL J n N M Re V w z

specific interfacial area (m−1 ) lift force coefficient (−) bubble diameter (m) Eötvös number (−) size distribution function (m−4 ) breakup and coalescence related variations of f (m−4 s−1 ) breakup, coalescence coefficients (−) lift force (kg m s−2 ) superficial velocity (m s−1 ) bubbles number density (m−3 ) number of velocity groups (−) number of sub-size groups (−) Reynolds number (−) velocity (m s−1 ) bubble velocity (m s−1 ) axial coordinate (m)

Subscripts B bubble l liquid g gas Greek symbols ˛ volumetric fraction (−)  density (kg m−3 )

(see Krepper et al., 2008). In case of vertical pipe flows, the boundary conditions are well defined and a cylindrical symmetry can be assumed. In addition, the bubbles move over a large distance forming in the result clear profiles. For this reason, pipe flow is a suitable subject for use in model development and validation. On the other hand, most of the flows relevant for practical applications are rather complex and show clear 3D effects. CFD codes are based on local models, which in principle should be independent of geometrical and scale effects. To demonstrate that the models proven for vertical pipe are also applicable for more complex flow situations, the present paper examines the simulation of experiments for the flow around the obstacle. 2. Models for polydispersed bubbly flows 2.1. General approach 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-specific continuity equations for volume fractions, and a separate set of momentum equations is solved for each phase. The exchange of momentum between phases is 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 wall lubrication and the turbulent dispersion force is described in detail in the next chapter. For the liquid a turbulence shear stress (SST) model according to Menter (1994) was applied. The SST model switches between the standard k-ε (for the flow away from walls) and the k-ω turbulence model (for the vicinity of walls) using a blending function, excluding the user influence modelling the near-wall conditions.

2373

Modelling the influence of the gas bubbles on the liquid turbulence the Sato’s eddy viscosity model for bubble induced turbulence was applied (Sato and Sekoguchi, 1975). The liquid viscosity  is modelled consistent of the laminar component L , the turbulence viscosity caused by the turbulence of the liquid phase T T = C L

k2 ε

with C = 0.09

(1)

and a bubble induced component: S = CS L ˛dB |V L − V G |

with CS = 0.6

(2)

2.2. Closure models for momentum exchange To model two-phase flow using the Euler/Euler approach, the momentum exchange between the phases must be considered. Apart from the drag that acts in flow direction, the so-called nondrag forces that act in the direction perpendicular to the flow direction must also be considered. Namely the lift force, the turbulence dispersion force, the virtual mass force, and the wall force play an important role in gas–liquid pipe flows. A more detailed discussion on force models and their influence on vertical pipe flow can be found at Lucas et al. (2007). The lift force FL considers the interaction of the bubble with the shear field of the liquid. For a single bubble, it reads 3

d  l ) × rot(w  l) g −w (w FL = −CL l 6

(3)

 g and where l is the liquid density, d the bubble diameter, and w  l are the bubble and liquid velocities, respectively. The classical w lift force, which has a positive coefficient CL , acts in the direction of decreasing liquid velocity. In the case of co-current upward pipe flow, this is the direction towards the pipe wall. Numerical (Ervin and Tryggvason, 1997) and experimental (Tomiyama et al., 1995) investigations have shown that the direction of the lift force changes its sign, if there is a substantial deformation of the bubble. Tomiyama (1998) investigated single bubble motion and derived the following correlation for the coefficient of the lift force from these experiments:



min[0.288 tan h(0.121Re), f (Eod )] Eod < 4 f (Eod ) for 4 < Eod < 10 −0.27 Eod > 10 with f (Eod ) = 0.00105Eo3d − 0.0159Eo2d − 0.0204Eod + 0.474 CL =

(4)

The coefficient CL depends on the modified Eötvös number Eod given by: Eod =

g(l − g )dh2 

(5)

Here dh is the maximum horizontal dimension of the bubble. It is calculated based on an empirical correlation for the aspect ratio developed by Wellek et al. (1966) with the following equation: dh = db

 3

1 + 0.163Eo0.757

(6)

and on the Reynolds number Re based on bubble diameter. Thus the sign of the lift coefficient, and consequently the direction of the lift force, depends on the bubble diameter. This behavior, originally found for single bubbles of air in glycerol was confirmed by a range of experiments also for gas–water polydispersed flows (e.g. Prasser et al., 2007), for steam–water flow and for different other fluid systems of gas–fluids too. For several flow configurations, this bubble size dependency of the lift force direction leads to the separation of small and large bubbles. This effect has been shown to be a key phenomenon for the development of the flow regime.

2374

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

The turbulent dispersion force is the result of the turbulent fluctuations of liquid velocity. It acts on the distribution of the gas volume fraction and can be evaluated from a single expression related to drag turbulent contribution. In previous investigations on vertical pipe flow, the model of Burns et al. (2004) was identified as most suitable: FTD = −

3CD t,l  w grad˛ 4db Pr l rel

(7)

To avoid the appearance of a maximum of the gas fraction at the wall Tomiyama et al. (1995) and Tomiyama (1998) proposed a wall force: d ˛ FW = −CW bubb 2



1 1 − y2 (D − y)2



2  nr l wrel

Fig. 1. Scheme of the standard MUSIG model: all size fractions representing different bubble sizes move with the same velocity field.

(8)

with the wall force coefficient



CW =

exp(−0.933Eo + 0.179)

1 ≤ Eo ≤ 5 (9)

for 0.007Eo + 0.04

5 ≤ Eo ≤ 33

This set of bubble forces was identified to reflect well the experimental results for fully developed upward vertical flow (Lucas et al., 2007; Krepper et al., 2008). For this reason, it was also applied in this work. The influence of the virtual mass force was found to be small. Thus, the virtual mass force was not considered in the simulations presented in this paper. 2.3. Bubble coalescence and breakup The bubble coalescence model takes into account the random collision processes between two bubbles. The applied model is based on the work of Prince and Blanch (1990). The bubble breakup model considers the collision between a liquid turbulent eddy of a certain characteristic size and a bubble (see Luo and Svendsen, 1996). The mechanistic models of coalescence and breakup are common to both methods, where the polydispersion in size of the bubble population is taken into account. We introduce corresponding “breakup” and “coalescence” coefficients FB and FC , as factors of the transition rates, which are used to scale the original correlations. 2.4. Population balance approach 2.4.1. The MUSIG model by Lo In principle, the Eulerian two-fluid approach can be extended to a multi-fluid approach, i.e. to simulate a continuous liquid phase and several gaseous dispersed phases solving the complete set of balance equations for each phase. However, investigations showed that for an adequate description of poly-dispersed flows, the inclusion of a population balance model requires decades of bubble size classes. In a CFD code, such a procedure is limited by the increased computational effort needed to obtain converged flow solutions. To solve this problem, the multiple size group model first implemented by the code developers in CFX-4 solves only one common momentum equation for all bubble size classes (homogeneous MUSIG model, see Lo, 1996, Fig. 1). Mathematically, the Multiple Size Group model (MUSIG) is based on the population balance method and the two-fluid modelling approach. The dispersed phase is divided into M size fractions. The population balance equation is used to describe the mass conservation of the size fractions and it accounts for the inter-fraction mass transfer caused by bubble coalescence and breakup. This approach allows a sufficient number of size fraction groups to be used for adequate coalescence and breakup calculations. A number of successful simulations of large-scale industrial multiphase flow problems have shown the applicability of the approach.

Fig. 2. Improvement of the MUSIG approach: the size fractions Mj are assigned to the velocity field Vj .

Nevertheless, the assumption also restricts its applicability to homogeneous dispersed flows, where the slip velocities of particles are almost independent of particle size and the particle relaxation time is sufficiently small with respect to inertial time scales. Thus, the asymptotic slip velocity must almost instantaneously be attained. Thus, the homogeneous MUSIG model described above fails to predict the correct phase distribution when heterogeneous particle motion becomes important. One example is the bubbly flow in vertical pipes where the non-drag forces play an essential role on the bubble motion. The lift force changes its sign when it is applied for large deformed bubbles, which are dominated by the asymmetrical wake (Tomiyama et al., 1995, see Section 2.1). For this reason, large bubbles tend to move to the pipe core region resulting in a core void maximum whereas a near-wall void peak is measured for small bubbles. The radial separation of small and large bubbles cannot be predicted by the homogeneous MUSIG model. This has been shown to be a key mechanism for the establishment of a certain flow regime. 2.4.2. New strategy: the Inhomogeneous MUSIG model A combination of the consideration of different dispersed phases and the algebraic multiple size group model was proposed to combine both the adequate number of bubble size classes for the simulation of coalescence and breakup and a limited number of dispersed gaseous phases to limit the computational effort (Krepper et al., 2005). The Inhomogeneous MUSIG model was developed in cooperation with ANSYS CFX and it has been available for use in CFX since its implementation in the version 10 release of the solver (Shi et al., 2004; Zwart et al., 2003; Frank et al., 2005 see Fig. 2). In the Inhomogeneous MUSIG model, the gaseous dispersed phase is divided into a number, N velocity groups (or phases), where each velocity group is characterized by its own velocity field. The overall bubble size distribution is further represented by dividing the bubble diameter range within each of the velocity groups j into a number, Mj (j = 1. . .N), bubble sub-size fractions. The population balance model considers bubble coalescence and bubble breakup, which are applied to the sub-size groups. Hence, the mass exchange between the sub-size groups can exceed the size ranges assigned

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

2375

Fig. 3. Sketch of the movable obstacle with driving mechanism—a half-moon shaped horizontal plate mounted on top of a toothed rod.

to the velocity groups, which results in the mass transfer terms between the different dispersed phase equations or velocity groups. The lower and upper boundaries of the intervals for the bubble size fractions can be controlled by either an equal distribution of bubble diameter, an equal bubble mass or it can be based on user definition of the bubble diameter ranges for each distinct bubble diameter fraction. The subdivision should be based on the physics of bubble motion for bubbles of different size, e.g. dissimilar behavior of distinctly sized bubbles with respect to lift force or turbulent dispersion. Extensive model validation calculations have shown that in most cases N = 2 or 3 velocity groups are sufficient in order to capture the main phenomena in bubbly or slug flows (Krepper et al., 2007, 2008). 3. The experiment—bubbly flow around an obstacle The experiment presented here was performed at the TOPFLOW facility of Forschungszentrum Dresden-Rossendorf (FZD). A large

pipe test section with a nominal diameter of DN200 was used to study the flow field around an asymmetric obstacle (see Fig. 3). This is a suitable test case for CFD code validation, since the obstacle creates a pronounced three-dimensional two-phase flow field. Curved streamlines, which form significant angles with the gravity vector, a recirculation zone in the wake and a flow separation at the edge of the obstacle are common in industrial components and installations. Wire-mesh technology (Prasser et al., 1998) is applied to measure the gas volume fraction and the gas velocity in different distances up- and downstream of the obstacle (Prasser et al., 2008). The sensor provides detailed data on the instantaneous flow structure with a high resolution in space and time. In particular, they allow the visualization of the structure of the gas–liquid interface (Prasser et al., 2002). For the obstacle experiment 3D fields of the distribution of the gas volume fraction around the obstacle and local bubble size distributions were obtained. These data are highly suitable data for CFD code validation. The 3D field of the vertical

Fig. 4. Comparison of time averaged contours calculated by CFX and measured values for the gas volume fraction (left) and the absolute value of the vertical component of the water velocity (right) up- and downstream of the obstacle in the air–water test run 096, JL = 1.017 m/s, JG = 0.0898 m/s (FB = FC = 0.05).

2376

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

component of the liquid velocity is also estimated by applying a special procedure to the data obtained from the wire-mesh sensor. This procedure as well as the experiment is described by Prasser et al. (2008). The tests are performed both for air/water and for steam/water. In the current paper, only selected adiabatic air/water tests were considered. 4. Numerical settings Pre-test calculations using ANSYS/CFX and applying a monodispersed bubble size approach were performed for the conditions of test run 074 (JL = 1.017 m/s, JG = 0.0368 m/s) (see Prasser et al., 2005; Frank et al., 2006, 2007). In the calculation, a fluid domain was modelled 1.5 m upstream and downstream the obstacle. Only half of the tube was simulated via the inclusion of a symmetry boundary condition set at the xz-plane of the geometry. In the present paper the Inhomogeneous MUSIG model approach was applied to air/water obstacle experiments run 096 (JL = 1.017 m/s, JG = 0.0898 m/s) and run 097 (JL = 1. 611 m/s, JG = 0.0898 m/s). The aim of the Inhomogeneous MUSIG approach is the simulation of the bubble size dependence of the momentum closure laws. The setup of the bubble size classes was therefore based on the lift force correlation developed by Tomiyama (see Eqs. (4–6)). To simulate the sign change of the lift coefficient the simulation distinguished two gaseous dispersed phases for small bubbles (CL > 0, db < 6 mm) and large bubbles (CL < 0, db > 6 mm). In the presented calculations for run 096 and run 097, 25 and 20 respectively subsize gas fractions representing equidistant bubble sizes up to 25 mm and 20 mm respectively were simulated, assigned to the two dispersed gaseous phases. The first six size groups were assigned to the first gaseous phase (or velocity group) and the remaining size groups were assigned to the second gaseous phase. The bubble size distribution measured at the largest upstream position was set as an inlet boundary condition for the calculation. 5. Comparison of measured and calculated results 5.1. The main phenomena observed The numerical results obtained by applying the Inhomogeneous MUSIG model have been compared to three-dimensional wiremesh sensor data in Fig. 4 (run 096). The water velocity and the total gaseous void fraction are represented. All qualitative details of the structure of the two-phase flow field around the obstacle could be reproduced. In the region directly behind the obstacle, a strong vortex of the liquid combined with the accumulation of gas is observed. The shape and extension of the recirculation zone for the calculated data agrees with the measured data. Upstream of the obstacle, a stagnation point with lower gas content is seen in experiment and calculation. Details, like the velocity and void fraction maxima above the gap between the circular edge of the obstacle and the inner wall of the pipe are also found in a good agreement between experiments and calculations. In the unobstructed cross-sectional part of the tube, a strong jet is established. Frank et al. (2007) discussed the details of the calculations. In the following section, the structure of the flow with respect to the population balance approach applied to the dispersed bubble phase is presented.

Fig. 5. Turbulence eddy dissipation (run 096) (CFX).

the applied bubble breakup model of Luo and Svendsen (1996), bubble breakup can be expected in regions showing high turbulent eddy dissipation. Fig. 5 presents maximum values of the numerically evaluated turbulent eddy dissipation rate at the edges of the obstacle. At the same time, the applied bubble coalescence model of Prince and Blanch (1990) indicates the strength of the coalescence in regions of bubble accumulation, i.e. in the wake behind the obstacle, as shown by the gas accumulation Fig. 4. Both bubble coalescence and bubble breakup are expected directly behind the obstacle (see distribution of turbulence dissipation Fig. 5). Thus, both effects might partially compensate each other. In Fig. 5 the position z = 0.08 m in the wake of the obstacle is marked. Fig. 6 shows the measured cross-sectionally averaged bubble size distributions upstream (z = −0.52 m), directly behind (z = 0.08 m) and downstream of the obstacle (z = 0.52 m). In the bubble accumulation zone at z = 0.08 m the cross-sectional average shows a shift towards larger bubbles. The calculated bubble size distributions (see Fig. 7 for run 096 and Fig. 8 for run 097), however, show a shift of the mean bubble diameter towards smaller bubbles in 10 cm behind the obstacle. In the calculations, the bubble breakup is apparently overestimated. For air/water flows in unobstructed vertical pipes, breakup or coalescence coefficients of 0.05

5.2. Phenomena in the wake of the obstacle 5.2.1. Bubble size distributions A more detailed understanding of the flow situation can be obtained by considering the bubble size distribution. According to

Fig. 6. Measured bubble size distribution for run 096.

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

2377

Fig. 7. Bubble size distributions for run 096 (JL = 1.017 m/s, JG = 0.0898 m/s) (CFX) (FB = FC = 0.05).

were found to yield good agreement to measurements. These coefficients were applied also here. 5.2.2. Bubbles streamlines The effects of the lateral motion of small and large bubbles can be revealed in more detail by studying bubble streamlines and by analyzing lift forces acting on bubbles of different sizes. On the one hand, the liquid flow carries the small bubbles into the region behind the obstacle (see Fig. 9). On the other hand, the air accumulation in the wake region leads to bubble coalescence and the generation of large bubbles as revealed by the analysis of experimental results. Fig. 10 illustrates lateral deviation due to lift force via the lift force vectors. Fig. 11 shows calculated and measured gas distributions up- and downstream of the obstacle resolved to bubble size classes for run 096. In calculations considering bubble fragmentation, this phenomenon is underestimated. In calculations considering bubble fragmentation, the effect of bubble coalescence is exceeded in this region by bubble fragmentation. Fig. 11a shows a solution for FB = 0, i.e. without the consideration of breakup.

Fig. 9. Streamlines for small (left) and large (right) bubbles (run 096) (CFX).

resulting phenomena become more pronounced with increasing water velocity. Therefore, run 097 is considered, where the liquid velocity was increased to JL = 1.611 m/s. Fig. 12 presents measured and calculated cross-sectional gas fraction distributions for this run. In the furthest downstream measured cross-section, an almost gas bubble free region is found in the experimental data. At the boundary of this almost gas free region an accumulation of bubbles is found (see Fig. 12, right side, z = 0.52 m). This effect is seen in almost all air/water measurements but not in the steam/water tests. The streamline representation of the calculations (Fig. 10 for run 096, which is quite similar to run 097), however, indicate large bubbles being directed into the jet caused by the lift force. This discrepancy between experiment and calculation can possibly be explained by the strong water velocity gradient near the jet.

5.3. Phenomena in the jet In the cross-sectional area beside the obstacle, a strong jet of high liquid velocity is established creating a strong shear flow. The

Fig. 8. Bubble size distributions for run 097 (JL = 1.611 m/s, JG = 0.0898 m/s) (CFX) (FB = FC = 0.05).

Fig. 10. Bubble lift force vectors for the different gas velocity groups (run 096) (CFX).

2378

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

Fig. 11. Calculated (left side) and measured (right side) gas distributions up- and downstream of the obstacle resolved to bubble size classes (run 096 JL = 1.017 m/s, JG = 0.0898 m/s, FB = 0, FC = 0.05).

Fig. 12. Calculated by CFX (left) and measured (right) gas cross fractional distributions downstream the obstacle (run 097 JL = 1.611 m/s, JG = 0.0898 m/s, FB = FC = 0.05) calculations (obstacle shown), distances at z = 0.08 m, 0.16 m, 0.25 m, 0.37 m and 0.52 m Measurements (obstacle in the upper left area) distances at z = 0.01 m, 0.015 m, 0.02 m, 0.04 m, 0.08 m, 0.16 m, 0.25 m and 0.52 m.

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

2379

Fig. 13. Calculated gas cross fractional distributions of small (left) and large (right) bubbles downstream the obstacle (run 097 JL = 1.611 m/s, JG = 0.0898 m/s, FB = FC = 0.05) calculations at distances from the obstacle of z = 0.08 m, 0.16 m, 0.25 m, 0.37 m and 0.52 m.

This strong shear flow induces bubble fragmentation, which is not yet considered in the model of Luo and Svendsen (1996). In the tests, the big bubbles migrated towards the jet, but they are fragmented at the relatively sharp boarder of this jet. Only a small fraction of the small bubbles created by this breakup process can enter the jet driven by the turbulent dispersion force. In the calculations, the furthest downstream plane from the obstacle shows an increase in the gas phase with larger bubbles and a decrease of the gas phase of smaller bubbles (see Fig. 13). The simulated total distribution (see Fig. 12) is therefore determined

Fig. 14. Modified coalescence coefficient (see Eq. (3)).

Fig. 15. Area of reduced bubble breakup.

2380

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

Fig. 16. Influence of the modified coalescence and reduced bubble breakup on the size decomposed gas volume fraction distribution.

by the simulated relation of the large and small bubbles fractions. However, the accumulation of small bubbles near the edge of the jet, which was found in the air/water tests, cannot be simulated via the described model approach applied here. 5.4. Improvement of the representation of bubble coalescence and breakup In Section 5.2, the sensitive overlay between bubble fragmentation and coalescence processes was shown. In this section, the effect of a recently published improvement to the coalescence model is investigated. In the literature, proposals for an increased coalescence coefficient for high void fractions can be found (e.g. Wang et al., 2005):



FC =



Min FCMAX

1 FC0 , FCMAX ˛0 − ˛





˛ ≤ ˛0

(10)

˛ ≥ ˛0

In the simulation, ˛0 was set to 0.225 and FC was limited to FCMAX = 100. This caused a significant increase in the bubble coalescence in the region directly behind the obstacle. Fig. 14 shows the distribution of the modified coalescence coefficient. At the same time, the turbulence eddy dissipation considered for the calculation of the bubble breakup was limited to occur at values less than 0.5 [m2 /s3 ] (see Fig. 15). These parameters were set to demonstrate the effect of the model modifications. Fig. 16 shows the effect on the distribution of the gas volume fraction for the conditions of run 096. In the original configuration of the Inhomogeneous MUSIG model where the breakup and coalescence coefficients are equal, the accumulation of small bubbles behind the obstacle is strongly overestimated (see Fig. 16a). By applying the described modifications to the coalescence coefficient and the limit of the eddy dissipation rate on breakup, a shift of the calculations in the right direction can be observed for the bubble size dependent gas distribution (see Fig. 11b for the measurements and Fig. 16b for the calculations). Fig. 17 presents the influence on the cross-sectional averaged bubble size distribution. Nevertheless, a satisfying agreement to the measurements cannot yet be achieved. 6. Summary and perspectives

Fig. 17. Influence of the modified coalescence and reduced bubble breakup on the cross-sectional averaged bubble size distribution.

In this study, the focus was directed on the model of the bubble size distribution in the numerical simulation of bubbly flows. A deeper understanding of the flow structure is possible when considering a more accurate characterization of the polydispersion. For upward two-phase flow in vertical pipes, the core peak in the crosssectional gas fraction distribution could be reproduced very well. For complex flows, the general three-dimensional structure of the flow could be well reproduced in the simulations. These test cases of pipe flow with an internal obstacle demonstrate the complicated relationship and interference between size dependent bubble migration, bubble coalescence and breakup effects in real flows. An advantage of the MUSIG model is that it does not impose any shape function for the bubble size distribution. The Inhomogeneous MUSIG model, using several velocity fields for the

E. Krepper et al. / Nuclear Engineering and Design 239 (2009) 2372–2381

bubbly phase, is able to deal with size separation of a locally polydispersed fluid in size population. These promising refinements of the two-fluid model for bubbly flows have shown their ability to recover consistent description of the lift force in the flow where complex flow structures are observed. While the closure models on bubble forces, which are responsible for the simulation of bubble migration, allow the explanation of experimentally observed hydrodynamic behavior of bubbles, clear deviations occur for bubble coalescence and breakup. In the simulations of the TOPFLOW experiments, the Luo and Svendsen model leads to an overestimation of the breakup that appears as negligible in the experiments. The presently applied models describing bubble breakup and coalescence could be proven as weak points in numerous CFD analyses. These bubble breakup and coalescence models depend to a large extent on the turbulence properties of the twophase flow, which were not measured and could not be validated in the pipe flow test cases. Therefore, further investigations are necessary to determine whether the multiphase flow turbulence models used here deliver appropriate and verifiable quantities that can be used for the description of bubble dynamics processes. Acknowledgements This study is carried out as part of current research projects funded by the German Federal Ministry of Economics and Labour, project numbers 150 1265 and 150 1329. The authors express their gratitude to the technical TOPFLOW (FZD) team. References Burns, A.D., Frank, T., Hamill, I., Shi, J.-M., 2004. The Favre averaged drag model for turbulent dispersion in eulerian multi-phase flows. In: 5th International Conference on Multiphase Flow, ICMF’04, Yokohama, Japan, May 30–June 4 (paper No. 392). Ervin, E.A., Tryggvason, G., 1997. The rise of bubbles in a vertical shear flow. J. Fluids Eng. 119, 443–449. Frank, T., Zwart, P.J., Shi, J.-M., Krepper, E., Rohde, U., 2005. Inhomogeneous MUSIG Model—a population balance approach for polydispersed bubbly flows. In: International Conference “Nuclear Energy for New Europe 2005, Bled, Slovenia, September 5–8, 2005. Frank, Th., Zwart, P.J., Krepper, E., Prasser, H.-M., Lucas, D., 2006. Validation of CFD models for mono- and polydisperse air–water two-phase flows in pipes. In: OECD/NEA International Workshop on The Benchmarking of CFD Codes for Application to Nuclear Reactor Safety (CFD4NRS), 05–09.09.2006, Garching, Deutschland. OECD/NEA, Garching, Germany.

2381

Frank, Th., Prasser, H.-M., Beyer, M., Al Issa, S., 2007. Gas–liquid flow around an obstacle in a vertical pipe—CFD simulation and comparison to experimental data. In: 6th Int. Conf. on Multiphase Flow Leipzig 2007, (paper 135). Krepper, E., Lucas, D., Prasser, H.-M., 2005. On the modelling of bubbly flow in vertical pipes. Nucl. Eng. Des. 235, 597–611. Krepper, E., Beyer, M., Frank, Th., Lucas, D., Prasser, H.-M., 2007. Application of a population balance approach for polydispersed bubbly flows. In: 6th Int. Conf. on Multiphase Flow Leipzig 2007, (paper 378). Krepper, E., Frank, Th., Lucas, D., Prasser, H.-M., Zwart, P.J., 2008. The Inhomogeneous MUSIG model for the simulation of polydispersed flows. Nucl. Eng. Des. 238, 1690–1702. Lo, S, 1996. Application of the MUSIG model to bubbly flows. In: AEAT-1096. AEA Technology. Lucas, D., Krepper, E., Prasser, H.-M., 2007. Use of models for lift, wall and turbulent dispersion forces acting on bubbles for poly-disperse flows. Chem. Sci. Eng. 62, 4146–4157. Luo, H., Svendsen, H.F., 1996. Theoretical model for drop and bubble break-up in turbulent flows. AIChE J. 42 (5), 1225–1233. Menter, F., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32 (8). Prasser, H.-M., Böttger, A., Zschau, J., 1998. A new electrode-mesh tomograph for gas–liquid flows. Flow Meas. Instrum. 9, 111–119. Prasser, H.-M., Krepper, E., Lucas, D., 2002. Evolution of the two-phase flow in a vertical tube—decomposition of gas fraction profiles according to bubble size classes using wire-mesh sensors. Int. J. Thermal Sci. 41, 17–28. Prasser, H.-M., Frank, T., Beyer, M., Carl, H., Pietruske, H., Schütz, P., 2005. Gas–liquid flow around an obstacle in a vertical pipe experiments and CFD simulation. In: Annual Meeting on Nuclear Technology, Nuremberg. Prasser, H.-M., Beyer, M., Carl, H., Gregor, S., Lucas, D., Pietruske, H., Schütz, P., Weiss, F.-P., 2007. Evolution of the structure of a gas–liquid two-phase flow in a large vertical pipe. Nucl. Eng. Des. 237, 1848–1861. Prasser, H.-M., Beyer, M., Al Issa, S., Carl, H., Pietruske, H., Schütz, P., 2008. Gas–liquid flow around an obstacle in a vertical pipe. Nucl. Eng. Des. 238, 1802–1819. Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and break-up in air-sparged bubble columns. AIChE J. 36 (10), 1485–1499. Sato, Y., Sekoguchi, K., 1975. Liquid velocity distribution in two phase bubble flow. Int. J. Multiphase Flow 2, 79–95. Shi, J.-M., Zwart, P.-J., Frank, T., Rohde U., Prasser, H.-M, 2004. Development of a multiple velocity multiple size group model for poly-dispersed multiphase flows. Annual Report of Institute of Safety Research. Forschungszentrum Rossendorf, Germany. Tomiyama, A., Sou, I., Zun, I., Kanami, N., Sakaguchi, T., 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: ICMF’98, 3rd Int. Conf. Multiphase Flow, Lyon, France, June 8–12, 1998, pp. 1–18. Wang, T., Wang, J., Yong, J., 2005. Theoretical prediction of flow regime transitions in bubble columns by the population balance model. Chem. Eng. Sci. 60, 6199– 6209. Wellek, R.M., Agrawal, A.K., Skelland, A.H.P., 1966. Shapes of liquid drops moving in liquid media. AIChE J. 12, 854–860. Zwart, P., Burns, A., Montavon, C., 2003. Multiple size group models. Technical Report. AEA Technology plc, November, 2003. CFX-5.7.