Science and Technology of Nuclear Installations
Volume 2009 (2009), Article ID 320738, 13 pages
doi:10.1155/2009/320738
Research Article

CFD Simulation of Polydispersed Bubbly Two-Phase Flow Around an Obstacle

E. Krepper,1 P. Ruyer,2 M. Beyer,1 D. Lucas,1 H.-M. Prasser,3 and N. Seiler2

1Institute of Safety Research, Forschungszentrum Dresden-Rossendorf (FZD), P.O. Box 510119, 01314 Dresden, Germany
2Institut de Radioprotection et de Sûreté Nucléaire, CE Cadarache, Bât. 700, BP 3, 13 115 Saint Paul lez Durance Cedex, France
3Institute of Energy Technology, ETH-Zürich, Sonneggstrasse 3, 8092 Zürich, Switzerland

Received 10 January 2008; Accepted 18 February 2008

Academic Editor: Iztok Tiselj

Copyright © 2009 E. Krepper et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

This paper concerns the model of a polydispersed bubble population in the frame of an ensemble averaged two-phase flow formulation. The ability of the moment density approach to represent bubble population size distribution within a multi-dimensional CFD code based on the two-fluid model is studied. Two different methods describing the polydispersion are presented: (i) a moment density method, developed at IRSN, to model the bubble size distribution function and (ii) a population balance method considering several different velocity fields of the gaseous phase. The first method is implemented in the Neptune_CFD code, whereas the second method is implemented in the CFD code ANSYS/CFX. Both methods consider coalescence and breakup phenomena and momentum interphase transfers related to drag and lift forces. Air-water bubbly flows in a vertical pipe with obstacle of the TOPFLOW experiments series performed at FZD are then used as simulations test cases. The numerical results, obtained with Neptune_CFD and with ANSYS/CFX, allow attesting the validity of the approaches. Perspectives concerning the improvement of the models, their validation, as well as the extension of their applicability range are discussed.

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 vapor of the liquid phase. The flow regimes found in vertical pipes are dependent on the void fraction of the gaseous phase, which varies, as void fraction increases, 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 monodispersed, an increase of the gas volume fraction leads to a broader bubble size distribution due to breakup and coalescence of bubbles. The exchange area for mass, momentum, or heat between continuous and dispersed phases thus cannot be simply modeled based on the knowledge of just the void fraction and a mean diameter. Moreover, the forces acting on the bubbles may depend on their individual size which is the case not only for drag but also for nondrag forces. Among the forces leading to lateral migration of the bubbles, that is, acting in normal direction with respect to the main drag force, 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 further coalescence of large bubbles migrating toward the pipe center into even larger Taylor bubbles or slugs.

An adequate modeling approach must consider all these phenomena. The paper presents two different approaches both being based on the Eulerian modeling framework. On one hand, a generalized inhomogeneous multiple size group (MUSIG) model was applied, for which the dispersed gaseous phase is divided into N inhomogeneous velocity groups (phases), each of these groups being subdivided into 𝑀 𝑗 bubble size classes. On the other hand, a moment density method is used to model the evolution of the bubble size distribution function along the flow, where main statistics of the distribution being described with a reduced set of transport equations. For both models, bubble breakup and coalescence processes are taken into account by appropriate models.

2. Model for the Local Size Distribution of a Bubble Population

2.1. Closure Models for Momentum Exchange and for Bubble Coalescence and Breakup

While modeling a two-phase flow using the Euler/Eulerian approach, the momentum exchange between the phases has to be considered. Apart from the drag acting in flow direction, the so-called nondrag forces acting mainly perpendicularly to the flow direction must be considered. Namely, the lift force, the turbulence dispersion force, the virtual mass force, and the wall force play an important role.

The turbulent dispersion force acts on smoothing the gas volume fraction distribution and can be evaluated either from a single expression related to drag turbulent contribution, [1] or deduced from a model for dispersed phase agitation contribution to the momentum balance, for example, the Tchen model [2]. To avoid the maximum gas fraction at the wall, Tomiyama et al. [3, 4] propose a wall force.

The lift force → 𝐹 𝐿 considers the interaction of the bubble with the shear field of the liquid. For a single bubble, it reads → 𝐹 𝐿 = − 𝐶 𝐿 𝜌 𝑙 𝜋 𝑑 3 6 ( → 𝑤 𝑔 − → 𝑤 𝑙 ) × r o t ( → 𝑤 𝑙 ) , ( 1 ) where 𝜌 𝑙 is the liquid density, 𝑑 the bubble diameter, and → 𝑤 𝑔 and → 𝑤 𝑙 are the bubble and liquid velocities, respectively. The classical lift force, which has a positive coefficient 𝐶 𝐿 , acts in the direction of decreasing liquid velocity. In case of cocurrent upwards pipe flow, this is the direction towards the pipe wall. Numerical [5] and experimental [3] investigations showed that the direction of the lift force changes its sign if a substantial deformation of the bubble occurs. Tomiyama [4] investigated single bubble motion and derived the following correlation for the coefficient of the lift force from these experiments: 𝐶 𝐿 = ⎧ ⎪ ⎨ ⎪ ⎩ m i n [ 0 . 2 8 8 t a n h ( 0 . 1 2 1 R e ) , 𝑓 ( E o 𝑑 ) ] , E o 𝑑 < 4 𝑓 ( E o 𝑑 ) , f o r 4 < E o 𝑑 < 1 0 − 0 . 2 7 E o 𝑑 > 1 0 w i t h 𝑓 ( E o 𝑑 ) = 0 . 0 0 1 0 5 E o 3 𝑑 − 0 . 0 1 5 9 E o 2 𝑑 − 0 . 0 2 0 4 E o 𝑑 + 0 . 4 7 4 . ( 2 ) This coefficient depends on the modified Eötvös number E o 𝑑 , which is formulated for the maximum horizontal bubble size in flow direction [6] and on the Reynolds number R e based on bubble size. The sign of the lift coefficient, and, consequently, the direction of the lift force, depends, thus, on the bubble size diameter. This behavior, originally found for single bubbles of air in glycerol, was also established by different experiments for gas-water polydispersed flows (e.g., [7]) and for different fluids too.

For several flow configurations, this bubble size dependency of the lift force direction can lead to the separation between small and large bubbles. This effect has been shown to be a key phenomenon for the development of the flow regime.

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 [8]. The bubble breakup model considers the collision between a liquid turbulent eddy of a certain characteristic size and a bubble (see [9]). Coalescence and breakup mechanistic models are common to both methods for taking into account the polydispersion in size of the bubble population. We introduce corresponding “breakup” and “coalescence” numerical coefficients 𝐹 𝐵 and 𝐹 𝐶 used to scale the original correlations.

2.2. Population Balance Approach
2.2.1. The MUSIG Model by Lo

In principle, the Eulerian two-fluid approach can be extended to simulate a continuous liquid phase and several gaseous dispersed phases solving the complete set of balance equations for each phase. The investigations, however, showed that for an adequate description of the gas volume fraction profile including a population balance model decades of bubble size classes would be necessary. In a CFD code, such a procedure is limited by the increased computational effort 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 [10], Figure 1). Mathematically, the multiple size group model (MUSIG) is based on the population balance method and the two-fluid modeling approach. The dispersed phase is divided into 𝑀 size fractions. The population balance equation is applied to describe the mass conservation of the size fractions taking into account the interfraction mass transfer resulting from bubble coalescence and breakup. This model approach allows a sufficient number of size fraction groups required for the coalescence and breakup calculation to be used and has found a number of successful applications to large-scale industrial multiphase flow problems.

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

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 can be considered to be attained almost instantaneously. 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 nondrag forces play an essential role on the bubble motion. The lift force was found to change its sign, when applied for large deformed bubbles, which are dominated by the asymmetrical wake ([3], see Section 2.1). The lift force in this case has a direction opposite to the shear-induced lift force on a small bubble. 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.2.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 [11]. The inhomogeneous MUSIG model was developed in cooperation with ANSYS CFX and is implemented in CFX since the version CFX-10 ([12–14] see Figure 2).

320738.fig.002
Figure 2: Improvement of the MUSIG approach: the size fractions 𝑀 𝑗 are assigned to the velocity field 𝑉 𝑗 .

In the inhomogeneous MUSIG model, the gaseous dispersed phase is divided into a number 𝑁 of the so-called velocity groups (or phases), where each of the velocity groups is characterized by its own velocity field. Further, the overall bubble size distribution is represented by dividing the bubble diameter range within each of the velocity groups 𝑗 in a number 𝑀 𝑗 , 𝑗 = 1 , … , 𝑁 , bubble subsize fractions. The population balance model, considering bubble coalescence or bubble breakup is applied to the subsize groups. Hence the mass exchange between the subsize groups can exceed the size ranges assigned to the velocity groups resulting in mass transfer terms between the different phases or velocity groups.

The lower and upper boundaries of bubble diameter intervals for the bubble size fractions can be controlled by either an equal bubble diameter distribution, an equal bubble mass distribution, or 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, for example, different behavior of differently sized bubbles with respect to lift force or turbulent dispersion. Extensive model validation calculations have shown that in most cases, 𝑁 = 2 or 3 , velocity groups are sufficient in order to capture the main phenomena in bubbly or slug flows [15, 16].

2.3. The Moment Density Method Approach

The moment density method proposes an alternative way to model polydispersion. This method allows to model the time and space evolution of a realistic distribution with the help of a very reduced number of transport equations, for example, Kamp et al. [17], Hill [18]. This is an interesting property, especially from a numerical point of view, with regard to alternative methods like population balance methods. In the present study, we restrict our attention to adiabatic bubbly flows of an incompressible gas inside a continuous liquid phase, and only consider a dispersion in bubbles size. We here present the basic formalism of the moment density method.

2.3.1. A Distribution Function Parameterized by Its Statistical Moments

The moment density method requires the use of an approximate representation of the bubble population thanks to a presumed shape continuous function for the bubble size distribution function 𝑓 , the one being thus totally determined by a finite number of parameters. The moment densities are intensive physical quantities characterizing mean properties of the population, like the density of bubble area or volume, the so-called volumetric interfacial area 𝑎 𝑖 , and the volumetric fraction 𝛼 of the dispersed phase. Bubbles being assumed as spherical and of diameter 𝑑 , it yields 𝑎 𝑖 = 𝜋 𝑑 2 𝑓 𝜕 𝑑 , 𝛼 = 𝜋 𝑑 3 6 𝑓 𝜕 𝑑 . ( 3 )

In this work, we introduce 𝑓 as being simply quadratic in the bubble size 𝑑 , ⎧ ⎪ ⎨ ⎪ ⎩ 𝑓 = 3 𝑛 ( → 𝑥 , 𝑡 ) 𝑑 4 𝑑 3 1 ( → 𝑥 , 𝑡 ) ( 2 𝑑 1 ( → 𝑥 , 𝑡 ) − 𝑑 ) i f 𝑑 ≤ 2 𝑑 1 , 0 e l s e w h e r e , ( 4 ) where 𝑛 is the local density number of bubbles ( m − 3 ) and 𝑑 1 is the mean diameter of the distribution. It is easy to show that 𝑛 , 𝑑 1 , 𝑎 𝑖 , and 𝛼 are related through: 2 𝑑 1 𝑎 𝑖 = 9 𝛼 , 2 4 3 𝜋 𝑛 𝛼 2 = 1 0 𝑎 3 𝑖 . ( 5 ) The knowledge of the two moment densities 𝑎 𝑖 and 𝛼 is thus sufficient to reconstruct the whole distribution function 𝑓 . The choice for these two parameters is dictated by their importance in the description of the two-phase flow, especially concerning the transfers between the phases.

2.3.2. Transport Equations for the Size Distribution

The transport of the distribution function 𝑓 obeys a Liouville-Boltzmann equation from which one can derive transport equations for the main moment densities as well as for their transport velocities [19]. The main idea of the method consists in solving a finite number of moment densities transport equations in the framework of the two-fluid model. This solving allows to determine the distribution parameters and thus the whole population evolution. Thus we consider the solving of the following transport equations: 𝜕 𝑎 𝑖 𝑎 𝜕 𝑡 + ∇ 𝑖 ⟨ → 𝑤 ⟩ 𝑎 𝑖 = 𝑑 2 𝜋 𝑑 𝑑 𝑑 𝑡 𝑓 𝜕 𝑑 + 𝜋 𝑑 2 ̇ 𝐹 𝑐 , 𝑏 𝜕 𝑑 , 𝜕 𝛼 𝜕 𝑡 + ∇ ⋅ 𝛼 ⟨ → 𝑤 ⟩ 𝛼 = 𝜋 𝑑 2 2 𝑑 𝑑 𝑑 𝑡 𝑓 𝜕 𝑑 , ( 6 ) where ̇ 𝐹 𝑐 , 𝑏 is the source term of the Liouville-Boltzmann equation. The transport velocities ⟨ → 𝑤 ⟩ 𝑎 𝑖 and ⟨ → 𝑤 ⟩ 𝛼 of 𝑎 𝑖 and 𝛼 , respectively, are assumed to be equal to the mean transport velocity of the bubbles → 𝑤 𝑔 (no correlation between dispersions in velocity and in size is taken into account at the current stage of development of the model). ⟨ → 𝑤 ⟩ 𝛼 is defined by 𝛼 ⟨ → 𝑤 ⟩ 𝛼 = 𝜋 𝑑 3 6 → 𝑤 𝑓 𝜕 𝑑 , ( 7 ) where → 𝑤 is the velocity of a bubble of size 𝑑 . It satisfies the transport equation: 𝜕 𝛼 ⟨ → 𝑤 ⟩ 𝛼 𝜕 𝑡 + ∇ ⋅ ( 𝛼 ⟨ → 𝑤 ⟩ 𝛼 ⟨ → 𝑤 ⟩ 𝛼 ) = − ∇ ⋅ ( 𝛼 ⟨ → 𝑤 → 𝑤 ⟩ 𝛼 ) + 𝜋 𝑑 3 6 𝑑 → 𝑤 𝑑 𝑡 𝑓 𝜕 𝑑 . ( 8 )

2.3.3. Closure Relations

The closure of the system of (3)–(8) concerns models (i) for coalescence and breakup events through ̇ 𝐹 𝑐 , 𝑏 , (ii) for the bubble size evolution along its trajectory 𝑑 𝑑 / 𝑑 𝑡 , (iii) for the hydrodynamic forces acting on individual bubbles, through 𝑑 → 𝑤 / 𝑑 𝑡 , as well as (iv) for the kinetic stress tensor ∇ ⋅ ( 𝛼 ⟨ → 𝑤 → 𝑤 ⟩ 𝛼 ) . The closure issue is based on the mechanistic model of all these phenomena at the scale of a single bubble. Source terms of the equations are then evaluated from the integral of individual contributions over the whole bubble population, whose size distribution function is known. The simplicity of the analytical expression (4) for 𝑓 allows to easily derive expressions for the integrals from the classical models described in Section 2.1. In this study, since the dispersed phase is considered as incompressible and incondensable, the bubble size evolution along its trajectory is zero. Detailed expressions for the source terms can be found in [20] .

It is worth pointing out that solving the system of (3)–(8), coupled with the system of balance equations for the continuous phase, is sufficient to describe the evolution of a broad spectrum of bubble sizes. This simplicity of the formulation is advantageous especially in terms of computational cost. Further extension of the model, taking into account the dispersion in velocity, can be considered using a similar formalism. For example, size-related velocity distribution can be introduced to model size-dependent bubble migration due to lift force, like it is done for the inhomogeneous MUSIG model. Moreover, this method has already been successful to consider more general particles velocity dispersion, for example, the work of Fox [21] concerning crossing trajectories.

3. The Experiment–Bubbly Flow Around an Obstacle

In the presented experiment performed at the TOPFLOW facility of Research center Dresden-Rossendorf (FZD), the large test section with a nominal diameter of DN200 was used to study the flow field around an asymmetric obstacle (see Figure 3). This is an ideal test case for CFD code validation, since the obstacle creates a pronounced three-dimensional two-phase flow field. Curved stream lines, 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.

320738.fig.003
Figure 3: Sketch of the movable obstacle with driving mechanism—a half-moon shaped horizontal plate mounted on top of a toothed rod.

The wire-mesh technology was applied to measure the gas volume fraction and the gas velocity in different distances up- and downstream the obstacle [22]. The sensor provides detailed data on the instantaneous flow structure with a high resolution in space and time. In particular, they allow visualizing the structure of the gas-liquid interface [23].

The tests were performed both for air/water and for steam/water. In the current paper, only adiabatic air/water tests were considered. The parameters are summarized in Table 1.

tab1
Table 1: Water and gas superficial velocities 𝐽 𝐿 and 𝐽 𝐺 of the tests investigated in this paper.

4. Numerical Settings

4.1. CFX and Population Balance Method

Pretest calculations using ANSYS/CFX and applying a monodispersed bubble size approach were performed for the conditions of test run 074 ( 𝐽 𝐿 = 1 . 0 1 7  m/s, 𝐽 𝐺 = 0 . 0 3 6 8  m/s) (see [24–26]). In the calculation, a fluid domain was modeled 1.5 m upstream and downstream the obstacle. Half of the tube including a symmetry boundary condition set at the 𝑥 𝑧 -plane of the geometry was simulated. In the present paper, the inhomogeneous MUSIG model approach was applied to air/water obstacle experiments run 096 ( 𝐽 𝐿 = 1 . 0 1 7  m/s, 𝐽 𝐺 = 0 . 0 8 9 8  m/s) and run 097 ( 𝐽 𝐿 = 1 . 6 1 1  m/s, 𝐽 𝐺 = 0 . 0 8 9 8  m/s). 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 2 dispersed gaseous phases. The first 6 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.

4.2. Neptune_CFD and Moment Density Method

Neptune_CFD code [27], which is based on the two-fluid model formulation, has been used to perform numerical simulations using the moment density method. The test run 074 has been simulated using a calculation domain corresponding to a 1 m long half pipe centered on the obstacle of around 150 000 cells that takes benefit of the symmetry of the experimental setup. Neptune_CFD calculations of test 074 are based on a 𝑘 - 𝜀 model for the liquid flow, the Tchen correlation being used to model the fluctuations related to bubbles turbulent motion. Bubble hydrodynamics model include both laminar and turbulent evaluations of drag, lift, and virtual mass forces, thanks to a drift model for the relative velocity. Inlet boundary conditions correspond to flat profiles for 𝛼 , 𝑎 𝑖 , and the velocities, allowing to recover the experimental surface averaged values for superficial velocities and mean bubble diameter. The surface averaged volumetric fraction and mean diameter are used to reconstruct the inlet bubble size distribution function. The numerical method is based on a pressure-based method. Mass, momentum, and energy equations are coupled by an iterative procedure within a time step. Spatial discretization is based on finite volume framework on unstructured meshes, all the variables being computed at the center of cells.

5. Comparison of Measured and Calculated Results

5.1. The Main Observed Phenomena

Both the steady-state ANSYS CFX calculations applying the inhomogeneous MUSIG model, and the Neptune_CFD calculations applying the moment density method, could reproduce all qualitative details of the flow structure of the two-phase flow field around the diaphragm. The structure of the flow for the here considered test cases 074, 096, and 097 are essentially similar. These different tests have been selected for the purpose of investigating certain phenomena which are more or less pronounced.

The numerical results have been compared to three-dimensional wire-mesh sensor data in Figure 4 (ANSYS CFX for the run 096) and Figure 5 (Neptune_CFD for the run 074). 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.

fig4
Figure 4: Comparison of time averaged values calculated by CFX (left) and measured (right) up- and downstream of the obstacle in the air-water test run 096, 𝐽 𝐿 = 1 . 0 1 7  m/s, 𝐽 𝐺 = 0 . 0 8 9 8  m/s, ( 𝐹 𝐵 = 𝐹 𝐶 = 0 . 0 5 ) .
320738.fig.005
Figure 5: Neptune_CFD numerical (left-hand side of each pair) versus TOPFLOW experimental (right-hand side of each pair) results for run 074. The two left-hand side pictures represent the dispersed phase volume fraction (scale from 0 to 6 ⋅ 1 0 − 2 ) across the symmetry plane of the pipe and the two right-hand side pictures the liquid velocity norm (scale from 0 to 2 ms−1).

Shortly, behind the obstacle a strong vortex of the liquid combined with the accumulation of gas is observed. The measured and calculated shape and extension of the recirculation area agree very well. Upstream 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. Main discrepancies between experimental and Neptune_CFD results concern the volumetric fraction upward (below) the obstacle and can be associated to the flat profiles used as inlet bottom boundary conditions.

The structure of the flow is studied in more detail in the following sections.

5.2. Phenomena in the Wake of the Obstacle
5.2.1. Size Distribution of the Bubble Population

More detailed understanding of the flow situation can be gained, considering the bubble size distribution. According to the applied bubble breakup model of Luo and Svendsen [9], bubble breakup can be expected in regions showing high-turbulent eddy dissipation. Figure 6 presents maximum values of the numerically evaluated turbulent eddy dissipation at the edges of the obstacle. At the same time, the applied bubble coalescence model of Prince and Blanch [8] indicates strong importance of coalescence in regions of bubble accumulation, that is, in the wake behind the obstacle. Both bubble coalescence (see gas accumulation shown in Figures 4 and 5) and bubble breakup (see distribution of turbulence dissipation Figure 6), which might partially compensate each other, are expected shortly behind the obstacle.

320738.fig.006
Figure 6: Turbulence eddy dissipation (run 096) (CFX).

Figure 7 presents measured cross-sectional averaged bubble size distributions upstream ( 𝑧 = − 0 . 5 2  m), shortly behind ( 𝑧 = 0 . 0 8  m) and downstream the obstacle ( 𝑧 = 0 . 5 2  m) for the run 096. The measurements show in the bubble accumulation zone at 𝑧 = 0 . 0 8  m the cross-sectional average shows a shift toward larger bubbles. The experimentally measured bubble size distribution of run 074 is essentially similar to the run 096. It can be seen on the volumetric fraction maps, Figure 8. Above the obstacle, where turbulent intensity is strong, breakup is not experimentally attested, contrarily to the Luo and Svedsen [9] model prediction, but coalescence occurs, leading to the formation of big bubbles (over 7 mm) that concentrate thereafter in the central part of the pipe.

320738.fig.007
Figure 7: Measured bubble size distribution for run 096.
fig8
Figure 8: TOPFLOW experimental volumetric fractions of air bubbles according to different size classes, extracted from [24].

Both ANSYS CFX and Neptune_CFD calculated bubble size distributions, however, show a shift of the mean bubble diameter toward smaller bubbles shortly behind the obstacle when both coalescence and breakup are taken into account. In the calculations, the bubble breakup is overestimated. The corresponding results are provided by Figure 9 for the run 096 and on Figure 10 for run 097. This disagreement was found not solvable by simple changes of breakup or coalescence coefficients, which were set here to 𝐹 𝐵 = 𝐹 𝐶 = 0 . 0 5 . Similar deviations would arise at other locations of the flow domain. Neptune_CFD results for the run 074 with standard coefficients 𝐹 𝐵 = 𝐹 𝐶 = 1 are shown in Figure 11 (the details can be found in [20]).

320738.fig.009
Figure 9: Bubble size distributions for run 096 ( 𝐽 𝐿 = 1 . 0 1 7  m/s, 𝐽 𝐺 = 0 . 0 8 9 8  m/s) (CFX) ( 𝐹 𝐵 = 𝐹 𝐶 = 0 . 0 5 ) .
320738.fig.010
Figure 10: Bubble size distributions for run 097 ( 𝐽 𝐿 = 1 . 6 1 1  m/s, 𝐽 𝐺 = 0 . 0 8 9 8  m/s) (CFX) ( 𝐹 𝐵 = 𝐹 𝐶 = 0 . 0 5 ) .
320738.fig.011
Figure 11: Neptune_CFD result for mean Sauter bubble diameter (m), both breakup and coalescence being taken into account ( 𝐹 𝐵 = 𝐹 𝐶 = 1 ) run 074.

This suggests to perform computations for which breakup is neglected, but that still consider the Prince and Blanch [8] coalescence model with standard coefficient. The corresponding Neptune_CFD results show that both bubble size repartition and order of magnitude are consistent with experimental results: big bubbles of more than 7 mm are created and accumulate in the recirculation zone just above the obstacle (see Figure 14). Big bubbles stay in the LHS of the pipe, the unobstructed cross-sectional part of the pipe is free of big bubbles, that is consistent with the experimental observation.

As partial conclusions, (i) the Luo and Svendsen breakup model tends to overestimate the breakup of TOPFLOW experimental tests, whereas (ii) the use of the Prince and Blanch model for coalescence within the framework of the moment density method allows to predict satisfying evolution of bubble size across the flow.

5.2.2. Bubbles Streamlines

More detailed effects of lateral motion of small and large bubbles can be revealed by studying bubble streamlines and by analyzing lift forces acting on bubbles of different size. On one hand, the liquid velocity flow carries the small bubbles into the region behind the obstacle (see Figures 12 and 14, right-hand side for the bubble streamlines). Lateral deviation due to lift force is illustrated on Figure 13 for the lift force arrows with the population balance method and Figure 14, right-hand side for the lift coefficient with the moment density method. 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. This phenomenon is underestimated in the calculations taking breakup into account.

320738.fig.012
Figure 12: Streamlines for small (left) and large (right) bubbles (run 096) (CFX).
320738.fig.013
Figure 13: Bubble lift force vectors for the different gas velocity groups (run 096) (CFX).
fig14
Figure 14: Neptune_CFD numerical results for local mean bubble diameter and corresponding streamlines colored by the equivalent mean lift coefficient (run 074) ( 𝐹 𝐵 = 0 ; 𝐹 𝐶 = 1 ) .

Caused by the lift force, large bubbles are redirected into the downstream jet (see Figure 13) once they can be formed in the wake by coalescence. The streamline representation (see Figure 12) clearly shows this phenomenon for large bubbles already present in the upstream flow.

In the actual version of the moment density method used in Neptune_CFD, dynamics of the bubble population is estimated using a single transport velocity. The averaged lift contribution takes into account both the bubble size distribution and the Tomiyama correlation. When a majority of bubbles are locally above the critical Eötvös number, the lift coefficient changes its sign. In this case, the direction of the lift force is changed, as it can be seen on Figure 14. As a consequence, the bubble streamlines above the obstacle deviate to the center of the pipe. This explains the low value of the dispersed phase volumetric fraction near the left-hand side part of the pipe wall above the obstacle for Neptune_CFD numerical results (see Figure 5). This result is fully consistent with experimental results.

As a partial conclusion concerning the bubble streamlines calculations, both methods showed their ability to consider the effect of bubble size on the lateral deviation of bubble streamlines due to lift force. This provides a more precise understanding and a more accurate prediction of bubbles repartition across the flow.

5.3. Phenomena in the Jet

In the cross-sectional area beside the obstacle, a strong jet is established creating strong shear flow. The resulting phenomena are more pronounced with increasing water velocity, like in run 097, where the liquid velocity was increased to 𝐽 𝐿 = 1 . 6 1 1  m/s. Figure 16 represents measured and ANSYS CFX calculated cross-sectional gas fraction distributions for this run. In the most downstream cross section of the measurements an almost gas bubble free region in the centre of the jet is found. Bubbles are collected at the edges of the jet in regions of largest water velocity gradient.

The streamline representation of the ANSYS CFX calculations, however, (Figure 12 for run 096, which is fairly similar to run 097) indicates large bubbles being directed into the jet caused by the lift force. This discrepancy between experiment and ANSYS CFX calculations can possibly be explained by the strong water velocity gradient near the jet. This strong shear flow induces bubble breakup which is not yet considered in the model of Luo and Svendsen [9]. In the tests, the big bubbles could migrate toward the jet, but be 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 by action of the turbulent dispersion force.

The gas distribution resolved by bubble size classes (see Figure 15) shows clearly the deviations between measurements and calculations. In the experiment, the gas accumulation behind the obstacle leads to a strong coalescence and the creation of large bubbles. In the CFX-calculations, this effect is underestimated.

fig15
Figure 15: Calculated by CFX (left) and measured (right) gas distributions up- and downstream of the obstacle resolved to bubble size classes (run 096 𝐽 𝐿 = 1 . 0 1 7  m/s, 𝐽 𝐺 = 0 . 0 8 9 8  m/s, 𝐹 𝐵 = 0 , 𝐹 𝐶 = 0 . 0 5 ).
320738.fig.016
Figure 16: Calculated by CFX (left) and measured (right) gas cross-fractional distributions downstream the obstacle (run 097 𝐽 𝐿 = 1 . 6 1 1  m/s, 𝐽 𝐺 = 0 . 0 8 9 8  m/s, 𝐹 𝐵 = 𝐹 𝐶 = 0 . 0 5 ). Calculations (obstacle shown), distances at 𝑧 = 0 . 0 8  m, 0.16 m, 0.25 m, 0.37 m, and 0.52 m. Measurements (obstacle in the upper left area). distances at 𝑧 = 0 . 0 1  m, 0.015 m, 0.02 m, 0.04 m, 0.08 m, 0.16 m, 0.25 m, and 0.52 m.

In Neptune_CFD calculations based on the moment density method the bubble breakup was neglected. As attested by Figure 14, left-hand side, in this case the bubbles present in the jet are small bubbles. Big bubbles created are concentrated in the center of the pipe where the largest value of volumetric fraction can be found (Figure 5). The structure of the gas repartition in the jet for both calculation and experiments is represented on Figure 17. In the upper part, the structure of the flow of the numerical calculation is consistent with experimental results. The most important deviation is observed in the region just above the obstacle ( 𝑧 = 0 . 0 8 m ) . A gain in accuracy of these results could be achieved by considering size-dependent velocity dispersion within the moment density method, in particular for bubble size-dependent lateral migration in the regions where there exists a strong mixing of different bubble sizes (mainly just above the obstacle as far as run 074 is concerned, see Figure 8).

320738.fig.017
Figure 17: Calculated by Neptune_CFD (left) and measured (right) gas volumetric fraction run 074 for different elevations (obstacle, symbolized in gray, is at 𝑧 = 0  m and covers the “upper” part of the half-pipe cross-section representation).

Nevertheless, the main trends of the bubble dynamics downward the obstacle are in agreement with experimental results (see Figure 16 for CFX and Figure 17 for Neptune_CFD).

6. Summary and Perspectives

In this study, we focused on the model of the bubble size distribution in the numerical simulation of bubbly flows. More deep 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 cross-sectional gas fraction distribution could be reproduced very well both by the moment density and by the MUSIG population balance methods. 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 internal obstacle demonstrate the complicated relationship and interference between size-dependent bubble migration, bubble coalescence, and breakup effects for real flows. With an appropriate given distribution function, the numerical effort of the moment density method is lower compared to the multiple bubble size group method (MUSIG). On the other hand, applying the MUSIG method the simulation of a flow situation allows to deal with more general shapes for the distribution function. Both methods enable to consider the effect of polydispersion in size on the bubble population dynamics, in particular on the evaluation of interphase momentum transfers associated with lift and drag. The inhomogeneous population balance model, using several velocity fields for the bubbly phase, is able to deal with size separation of a locally polydispersed in size population, whereas the moment density method accounts for local diversity in bubble hydrodynamics thanks to a single-averaged contribution of interphase transfers. These two promising refinements of the two-fluid model for bubbly flows have shown their ability to recover consistent description of the lift force in the upper part of 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 to explain the bubbles hydrodynamic behavior observed experimentally, clear deviations occur for bubble coalescence and breakup. Both methods based on similar coalescence and breakup models lead to the same conclusion: 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. On the other hand, the coalescence model of Prince and Blanch seems able to recover the correct bubble size, if used solely, as attested by corresponding Neptune_CFD calculations. For both methods, 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 two-phase 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 currently used multiphase flow turbulence models deliver appropriate and verifiable quantities that can be used for the description of bubble dynamics processes.

Extensions of both the moment density method and the MUSIG method to nucleate boiling regime numerical simulation are in progress. This includes the phenomena of compressibility, phase-change, and wall nucleation. To model the bubble size-dependent lateral migration phenomenon, the moment density method should also include a model for the bubble velocity distribution. This can be done using a similar formalism to the present model for bubble size distribution function.

Nomenclature
𝑎 𝑖 :Specific interfacial area [m−1]
𝐶 𝐿 :Lift force coefficient [-]
𝑑 :Bubble diameter [m]
𝑑 :Lagrangian derivative
E o : Eötvös number [-]
𝑓 :Size distribution function [m−4]
̇ 𝐹 𝑐 , 𝑏 :Breakup and coalescence related variations of 𝑓 [m−4s−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 [-]
R e :Reynolds number [-]
𝑉 :Velocity [m s−1]
𝑤 :Bubble velocity [m s−1]
𝑧 :Axial coordinate [m]
𝑙 :Liquid
𝑔 :Gas
𝛼 :Volumetric fraction [-]
𝜌 : Density [kg m−3].

Acknowledgments

Part of this study is carried out at the Institute of Safety Research of the FZD as a part of current research projects funded by the German Federal Ministry of Economics and Labour, Project nos. 150 1265 and 150 1271. The other part of this study carried out at the Institut de Radioprotection et de Sûreté Nucléaire (IRSN) has been achieved in the framework of the Neptune project, financially supported by CEA (Commissariat à l'Énergie Atomique), EDF (Électricité de France), IRSN, and AREVA-NP. The authors express their gratitude to the technical TOPFLOW (FZD) team.

References

  1. A. D. Burns, T. Frank, I. Hamill, and J.-M. Shi, “The favre averaged drag model for turbulent dispersion in Eulerian multi-phase flows,” in Proceedings of the 5th International Conference on Multiphase Flow (ICMF '04), pp. 1–17, Yokohama, Japan, May-June 2004, paper no. 392.
  2. E. Deutsch and O. Simonin, “Large eddy simulation applied to the motion of particles in stationary homogeneous fluid turbulence,” in Turbulence Modification in Multiphase Flows, vol. 110, pp. 35–42, ASME FED, Portland, Ore, USA, 1991.
  3. E. A. Ervin and G. Tryggvason, “The rise of bubbles in a vertical shear flow,” Journal of Fluids Engineering, vol. 119, no. 2, pp. 443–449, 1997.
  4. R. O. Fox, “Introduction and fundamentals of modelling approaches for polydisperse multiphase flows,” in Multiphase Reacting Flows: Modelling and Simulation, D. L. Marchisio and R. O. Fox, Eds., vol. 492 of CISM Courses and Lectures, pp. 1–40, Springer, New York, NY, USA, 2007.
  5. T. Frank, P. J. Zwart, J.-M. Shi, E. Krepper, D. Lucas, and U. Rohde, “Inhomogeneous MUSIG model—a population balance approach for polydispersed bubbly flows,” in Proceedings of the International Conference on Nuclear Energy for New Europe, Bled, Slovenia, September 2005.
  6. Th. Frank, P. J. Zwart, E. Krepper, H.-M. Prasser, and D. Lucas, “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 '06), Garching, Germany, September 2006, paper no. B6-32.
  7. Th. Frank, H.-M. Prasser, M. Beyer, and S. Al Issa, “Gas-liquid flow around an obstacle in a vertical pipe—CFD simulation and comparison to experimental data,” in Proceedings of the 6th International Conference on Multiphase Flow (ICMF '06), Leipzig, Germany, July 2007, paper no. 135.
  8. A. Guelfi, D. Bestion, M. Boucker, et al., “NEPTUNE: a new software platform for advanced nuclear thermal hydraulics,” Nuclear Science and Engineering, vol. 156, no. 3, pp. 281–324, 2007.
  9. D. P. Hill, The computer simulation of dispersed two-phase flows, Ph.D. thesis, Imperial College of Science, Technology and Medecine, University of London, London, UK, 1998.
  10. A. M. Kamp, A. K. Chesters, C. Colin, and J. Fabre, “Bubble coalescence in turbulent flows: a mechanistic model for turbulence-induced coalescene applied to microgravity bubbly pipe flow,” International Journal of Multiphase Flow, vol. 27, no. 8, pp. 1363–1396, 2001.
  11. E. Krepper, Th. Frank, D. Lucas, H.-M. Prasser, and P. J. Zwart, “Inhomogeneous MUSIG model-a population balance approach for polydispersed bubbly flows,” in Proceedings of the 6th International Conference on Multiphase Flow (ICMF '06), Leipzig, Germany, July 2007, paper no. 375.
  12. E. Krepper, M. Beyer, Th. Frank, D. Lucas, and H.-M. Prasser, “Application of a population balance approach for polydispersed bubbly flows,” in Proceedings of the 6th International Conference on Multiphase Flow (ICMF '06), Leipzig, Germany, July 2007, paper no. 378.
  13. E. Krepper, D. Lucas, and H.-M. Prasser, “On the modelling of bubbly flow in vertical pipes,” Nuclear Engineering and Design, vol. 235, no. 5, pp. 597–611, 2005.
  14. S. Lo, June 1996, Application of the MUSIG model to bubbly flows, AEAT-1096, AEA Technology.
  15. H. Luo and H. F. Svendsen, “Theoretical model for drop and bubble break-up in turbulent flows,” AIChE Journal, vol. 42, no. 5, pp. 1225–1233, 1996.
  16. H.-M. Prasser, A. Böttger, and J. Zschau, “A new electrode-mesh tomograph for gas-liquid flows,” Flow Measurement and Instrumentation, vol. 9, no. 2, pp. 111–119, 1998.
  17. H.-M. Prasser, E. Krepper, and D. Lucas, “Evolution of the two-phase flow in a vertical tube—decomposition of gas fraction profiles according to bubble size classes using wire-mesh sensors,” International Journal of Thermal Sciences, vol. 41, no. 1, pp. 17–28, 2002.
  18. H.-M. Prasser, T. Frank, M. Beyer, H. Carl, H. Pietruske, and P. Schütz, “Gas-liquid flow around an obstacle in a vertical pipe experiments and CFD simulation,” in Proceedings of the Annual Meeting on Nuclear Technology, vol. 16, Aachen, Germany, May 2005.
  19. H.-M. Prasser, M. Beyer, H. Carl, et al., “Evolution of the structure of a gas-liquid two-phase flow in a large vertical pipe,” Nuclear Engineering and Design, vol. 237, no. 15–17, pp. 1848–1861, 2007.
  20. M. J. Prince and H. W. Blanch, “Bubble coalescence and break-up in air-sparged bubble columns,” AIChE Journal, vol. 36, no. 10, pp. 1485–1499, 1990.
  21. P. Ruyer, N. Seiler, M. Beyer, and F. P. Weiss, “A bubble size distribution model for the numerical simulation of bubbly flows,” in Proceedings of the 6th International Conference on Multiphase Flow (ICMF '06), Leipzig, Germany, July 2007, paper no. 484.
  22. J.-M. Shi, P.-J. Zwart, Th. Frank, U. Rohde, and H.-M. Prasser, “Development of a multiple velocity multiple size group model for poly-dispersed multiphase flows,” 2004, Annual Report of Institute of Safety Research. Forschungszentrum Rossendorf, Dresden, Germany.
  23. O. Simonin, “Continuum modelling of dispersed two-phase flows,” in Combustion and Turbulence in Two Phase Flows, VKI Lecture Series 1996-02, Von Karman Institute for Fluid Dynamics, Brussels, Belgium, 1996.
  24. A. Tomiyama, I. Sou, I. Zun, N. Kanami, and T. Sakaguchi, “Effects of Eötvös number and dimensionless liquid volumetric flux on lateral motion of a bubble in a laminar duct flow,” in Advances in Multiphase Flow, pp. 3–15, Elsevier Science, Amersterdam, The Netherlands, 1995.
  25. A. Tomiyama, “Struggle with computational bubble dynamics,” in Proceedings of the 3rd International Conference on Multiphase Flow (ICMF '98), pp. 1–18, Lyon, France, June1998.
  26. R. M. Wellek, A. K. Agrawal, and A. H. P. Skelland, “Shapes of liquid drops moving in liquid media,” AIChE Journal, vol. 12, no. 5, pp. 854–862, 1966.
  27. P. A. Zwart, A. Burns, and C. Montavon, “Multiple size group models,” AEA Technology plc, Oxfordshire, UK, November 2003.