Population balance modeling for air–water bubbly flow in a vertical U-bend

Abstract

Bubbly flow in U-bend is widely encountered in two-phase flow systems because of its compactness and high heat transfer coefficient. The modeling of phase distributions, velocity fields, and interfacial area concentration in the U-bend is crucial for the analysis of mass, momentum, and energy transportation processes in the equipment. However, this subject has not received enough attention yet. In this paper, the combination of population balance model and two-fluid model was used in the simulation of air–water bubbly flow in a U-bend with 24 mm inner diameter and 96 mm curvature. The homogeneous multiple size group model was used to solve the population balance equation and reconstruct the bubble size distribution function. The phase distribution at 0°, 90°, and 180° was predicted and the results showed that the superficial velocities of gas and liquid phase were the control parameters. Under higher gas superficial velocity, the buoyant force is dominant and makes the bubbles concentrate on the outer side of the tube wall; while under lower gas superficial velocity, the centrifugal force is dominant and makes the bubbles concentrate on the inner side of the tube wall. These results met well with the experimental results of Usui.

Introduction

The vertical U-bend is widely used as the basic heat transfer elements of heat exchanger, boiler, steam generator, and evaporator, in the plant of chemical, energy, and nuclear industries, because of its compactness and high efficiency. Frequently, to achieve high heat transfer coefficient, the equipment are designed to work in subcooled boiling region and the flow behaves as a bubbly flow. Under such flow condition the behavior of bubbles influences the local heat and mass transfer processes significantly. It is necessary to conduct accurate predictions of the flow characteristics of bubbly flow in the vertical U-tube to improve the economic efficiency and safety of the plant.

The numerical studies of bubbly flow were mainly concentrated on vertical or horizontal straight tubes.^1,2 In straight tubes the movement of bubbles was primarily influenced by the gravity force, which controls the gas–liquid relative velocity (for vertical tube) or phase distribution (for horizontal tube). Nevertheless, in the vertical U-bend the centrifugal force is exerted on the fluids in the direction perpendicular to the velocity and influences the phase distribution dramatically. Under the coaction of gravity and centrifugal force, the flow phenomenon becomes very complicated, which increases the difficulty in modeling bubbly flow in vertical U-bend in comparison with straight flow.

Several experimental studies were performed on the gas–liquid bubbly flow in the U-bend. Usui et al.^3,4 investigated the transition of flow patterns in an invert U-bend with tube diameter of 16 and 24 mm and curvature of the bend of 90–180 mm. They also measured the phase distributions of different cross section in bubbly flow and annular flow region using an electrical conductivity probe. Based on the experimental results the authors discussed the influence of gravity and centrifugal force on the phase distribution. Hoang and Davis ⁵ measured the void fraction, local gas velocity, and bubble size along an inverted U-bend with tube diameter of 50.8 mm and bend curvature of 50.8 and 76.2 mm. Arlindo and Franca ⁶ performed both experimental and numerical analysis on the air–water bubbly flow in a U-bend with square cross section. The phenomena like the wall peaking, happening in vertical upward flows; the core peaking, taking place in downward flows; and the action of the centrifugal force, setting the phase distribution along the bend were measured. Meng et al. ⁷ visually observed the flow boiling and two-phase flow in U-bend of a vertical downward hairpin tube with inner diameter of 6 mm. The bubble dynamical behavior was divided into four zones based on flow pattern and bubble size. The phase distribution in the U-bend was not tested directly in the experiments.

Few numerical simulations were carried out for the bubbly flow inside a U-bend. Wu et al. ⁸ simulated refrigerant flow boiling in a serpentine round tube using Eulerian two-fluid model and a phase-change model for the mass transfer. In the simulation the orientation of the U-bend was horizontal and the phase distribution was not confirmed by the test. Arlindo and Franca ⁶ performed numerical simulation using two-fluid model and interphase momentum transfer correlations. The numerical results of phase distribution met well with the experiments but the bubble diameter was an input to the numerical model based on the test results instead of computational results.

In recent years, in order to predict the bubble size distribution (BSD) and avoid arbitrary mean bubble diameter input, the population balance method was introduced and applied to simulation of bubbly flow in various channels. In this approach the bubble’s characteristics, such as size, surface area, composition, and so forth, are treated as “internal” coordinate. The bubble numbers of each coordinate were controlled by the bubble birth or death mechanisms caused by bubble coalescence and breakage, and described by the population balance equation (PBE). By solving PBE the BSD can be traced automatically. Ekambara et al., ⁹ Li et al., ¹⁰ Yeoh et al., ¹ and Yeoh and Tu ² have adopted this approach to investigate the internal phase distribution of bubbly flow in vertical or horizontal straight tubes.

In this paper we focus on the prediction of the phase distribution and phase velocities in a vertical U-bend using the population method and two-fluid model. The homogeneous multiple size group (MUSIG) model, introduced by Lo and Zhang, ¹¹ was used to solve the PBE and reconstruct the BSD function. The continuous BSD function was approximated by a set of discrete size fractions; the mass conservation of each size fraction is balanced by the interfraction mass transfer due to the mechanisms of bubble–bubble interactions due to coalescence and breakage processes. The phase distributions at different positions along the bend were calculated; the numerical results met well with the test results of Usui et al. ³ This means the population method can offer a powerful approach to describe the bubbly flow in a U-bend. Based on simulation results of phase distribution characteristics under different superficial velocities, the influences of centrifugal force and gravity force were discussed.

Numerical models

Two-fluid model

The two-fluid model considers gas phase and liquid phase separately and is expressed in terms of two sets of conservation equations governing the balance of mass, momentum, and energy in each phase. For isothermal two-phase flow, the energy equation should be eliminated, hence the control equations can be written as:

Continuity equation of liquid phase

$$\frac{\partial }{\partial t}\left({\rho }_l{\alpha }_l\right)+\nabla \cdot \left({\rho }_l{\alpha }_l{\mathbf{u}}_l\right)=S_l$$

(1)

Continuity equation of gas phase

$$\frac{\partial }{\partial t}\left({\rho }_g{\alpha }_g\right)+\nabla \cdot \left({\rho }_g{\alpha }_g{\mathbf{u}}_g\right)=S_g$$

(2)

Momentum equation of liquid phase

$$\begin{array}{c}\frac{\partial }{\partial t}\left({\rho }_l{\alpha }_l{\mathbf{u}}_l\right)+\nabla \cdot \left({\rho }_l{\alpha }_l{\mathbf{u}}_l{\mathbf{u}}_l\right)=-{\alpha }_l\nabla p+p_l\mathbf{g}\\ +\nabla \cdot \left[{\alpha }_l{\mu }_{e,l}\left(\nabla {\mathbf{u}}_l+{\nabla }^T{\mathbf{u}}_l\right)\right]+F_{gl}\end{array}$$

(3)

Momentum equation of gas phase

$$\begin{array}{c}\frac{\partial }{\partial t}\left({\rho }_g{\alpha }_g{\mathbf{u}}_g\right)+\nabla \cdot \left({\rho }_g{\alpha }_g{\mathbf{u}}_g{\mathbf{u}}_g\right)=-{\alpha }_g\nabla p+p_g\mathbf{g}\\ +\nabla \cdot \left[{\alpha }_g{\mu }_{e,g}\left(\nabla {\mathbf{u}}_g+{\nabla }^T{\mathbf{u}}_g\right)\right]+F_{\lg }\end{array}$$

(4)

where the subscripts g and l represent the parameters of gas phase and liquid phase, respectively; $F_{lg}$ and $F_{gl}$ are the interface forces between gas and liquid phase ( $F_{lg}=F_{gl}$ ) and can be expressed as

$$F_{lg}=F_{lg}^D+F_{lg}^L+F_{lg}^{WL}+F_{lg}^{TD}$$

(5)

where $F_{lg}^D$ is the drag force, $F_{lg}^L$ is the lift force, $F_{lg}^{WL}$ is the wall lubrication force, and $F_{lg}^{TD}$ is the turbulence dispersion force.

The drag force is induced by the relative movement between the phases and has the same direction with relative velocity, which can be obtained by Ishii and Zuber model

$$F_{\lg }^D=\frac{1}{8}{C}_D{a}_{if}{\rho }_l\left|{\mathbf{u}}_g-{\mathbf{u}}_l\right|({\mathbf{u}}_g-{\mathbf{u}}_l)$$

(6)

where $C_D$ is the drag force coefficient and $a_{if}$ is the interfacial area concentration. $C_D$ has been correlated for different Reynolds number regions for individual bubbles by Ishii and Zuber.

Lift force is relative to the vorticity field and directed perpendicular to the relative velocity of the phases. In this paper the lift force was calculated by Tomiyama model

$$F_{lg}^L=-C_L{\alpha }_g{\rho }_l\left({\mathbf{u}}_g-{\mathbf{\text{u}}}_l\right)\times \nabla \times {\mathbf{\text{u}}}_l$$

(7)

where $C_L$ is the lift force coefficient, determined by the Reynolds number, Eotvos number, and bubble size.

Wall lubrication force reflects the influence of surface tension of the tube wall, which exerts on the bubbles in the direction pointing to the tube axis. The correlation proposed by Frank was adopted in this paper

$$F_{lg}^{WL}=-C_{WL}{\alpha }_g{\rho }_l{\left|{\mathbf{\text{u}}}_g-{\mathbf{\text{u}}}_l\right|}^2{\mathbf{\text{n}}}_W$$

(8)

where $C_{WL}$ is given by

$$C_{WL}=C_W\left(Eo\right)\cdot \max \left\{0,\frac{1}{C_{WD}}\cdot \frac{1-\frac{y_w}{C_{WC}{d}_p}}{y_W\cdot {\left(\frac{y_W}{C_{WC}{d}_p}\right)}^{p-1}}\right\}$$

(9)

where $C_{WC}$ , $C_{WD}$ , and $p$ are the coefficients of this model and their values vary for different problems. The Antal model (with calibration factor C_W1 and C_W2) and Tomiyama model were also tested in comparison with Frank model, whose details can be found in Team ANSYS-CFXD. ¹²

Turbulent dispersion force is induced by turbulent flow and has the effect to decrease the void fraction gradient. The correlation proposed by Lopez de Bertodano was adopted

$$F_{lg}^{TD}=-C_{TD}{\rho }_l\nabla {\alpha }_l{k}_l$$

(10)

where $C_{TD}$ is turbulence dispersion coefficient, whose value is between 0.1 and 0.5.

Population balance method

For bubbly flow, the mass, momentum, and energy transportations between the phases were heavily affected by the area concentrations or bubble size, as shown by equation (10). Therefore, it is crucial to predict the local BSD accurately for a successful modeling of bubbly flow. This goal can be accomplished by the population balance method through solving the control equation of particle size distribution (PSD) function. In this paper the MUSIG model was adopted. As one kind of population balance method, MUSIG model approximates the continuous PSD with a series number of discrete bubble size classes. For each class a scalar PSD equation is solved to obtain the population changes caused by bubble coalescence and breakage. For ith class the bubble’s PSD is described by

$$\frac{\partial n_i}{\partial t}+\nabla \cdot \left({\mathbf{\text{u}}}_g{n}_i\right)=B_C+B_B-D_C-D_B$$

(11)

where $n_i$ is the bubble number density; $B_C$ and $B_B$ are the bubble number density increase rate induced by coalescence and breakage, respectively; $D_C$ and $D_B$ are the bubble number density decrease rate induced by coalescence and breakage, respectively. The following relationships were provided

$$B_{C,i}=\frac{1}{2}{\displaystyle \sum _{j=1}^i\left(a_{i-j,j}{n}_{i-j}{n}_j\right)}$$

(12)

$$D_{C,i}=n_i{\displaystyle \sum _{j=1}^{N-i}\left(a_{i,j}{n}_j\right)}$$

(13)

$$B_{B,i}=n_i{\displaystyle \sum _{j=i+1}^N\left(g_{j,i}{n}_j\right)}$$

(14)

$$D_{B,i}=n_i{\displaystyle \sum _{j=1}^{i-1}\left(g_{i,j}\right)}$$

(15)

where $a_{i,j}$ is the coalescence rate between the bubbles of size class i and j and $g_{j,i}$ the breakage rate of bubbles from class j to class i. In this paper the coalescence rate was calculated by the correlation from Prince and Blanch, ¹³ and the breakage rate was calculated by the correlation from Luo and Svendsen, ¹⁴ respectively

$$a_{i,j}=F_C\frac{\pi }{4}{\left[d_i+d_j\right]}^2{\left(u_{ti}^2-u_{tj}^2\right)}^{0.5}\exp \left(-\frac{t_{ij}}{{\tau }_{ij}}\right)$$

(16)

$$\begin{array}{c}{g}_{i,j}=F_{B}C\left(1-{\alpha }^g\right){\left(\frac{{\epsilon }^l}{d_j^2}\right)}^{\frac{1}{3}}\\ {\displaystyle {\int }_{{\xi }_{\min }}\frac{{\left(\xi +1\right)}^2}{{\xi }^{\frac{11}{3}}}\exp \left(-\frac{12c_f\sigma }{\beta {\rho }^l{\left({\epsilon }^l\right)}^{\frac{2}{3}}{d}^{\frac{5}{3}}{\xi }^{\frac{11}{3}}}\right)\text{d}\xi }\end{array}$$

(17)

where $F_C$ is the coalescence calibration factor, $F_B$ is the breakup calibration factor, ${\tau }_{ij}$ is the contact time for two bubbles, and $t_{ij}$ is the time required for two bubbles to coalesce; $\xi =\frac{\lambda }{d}$ is the size ratio between an eddy and a particle; $C$ and $\beta$ are determined from fundamental consideration of bubbles breakup in turbulent dispersion systems to be 0.923 and 2.0; $c_f=f_{BV}^{\hspace{0.17em}2/3}+{\left(1-f_{BV}\right)}^{2/3}-1$ , $f_{BV}$ is stochastic breakup volume fraction of the daughter size distribution.

Numerical details

The parameters of the inverted U-bend analyzed in this paper are illustrated in Figure 1. The U-bend was invertedly placed and the mixture of air and water flows upward toward the bend, then, after passing the bend, flows downward toward the outlet. The inside diameter of the tube is 24 mm, while the radius of curvature is 96 mm. The lengths of the straight part at the inlet and outlet of the bend were set to be 1 and 0.3 m, respectively, to get fully developed flow status and prevent reverse flow. The working pressure and temperature of the fluids are 1 atm and 25°C, respectively. Fully dispersed bubbly flow only occurs under the conditions with large water superficial velocity and small air superficial velocity, according to the result of Usui et al.³ Therefore, the following cases were investigated in this paper (Table 1).

Figure 1.

Structure of the inverted U-bend.

Table 1.

Parameters of test section.

Case number	Superficial velocity of water (m/s)	Superficial velocity of air (m/s)
1	1	0.131
2	4	0.131
3	4	0.5
4	4	1.0

The mesh was generated by ICEM software, which has a total of 483,000 nodes, as shown in Figure 2. We have carried out mesh irrelevant study and found that a finer mesh cannot improve the calculation obviously and this mesh was used throughout the paper. Commercial CFD software CFX14 was used to perform numerical simulation. The k–e model and scalable wall function was adopted for the turbulence modeling. For two-phase flow simulation, the homogenous MUSIG model was adopted and treated the air phase as polydispersed fluid. The bubble size was dived into 10 classes, as shown in Table 2. The inlet boundary condition was selected as constant velocity condition for both water and air phase, while the outlet boundary was selected as constant pressure condition (zero relative pressure).

Figure 2.

Mesh of the U-bend.

Table 2.

Bubble diameters of each bubble class.

Bubble classes	Diameters (mm)
1	0.69
2	1.67
3	2.65
4	3.63
5	4.61
6	5.59
7	6.57
8	7.55
9	8.53
10	9.51

Results and discussion

Sensitivity analysis of interphase momentum transfer model

As mentioned in “Numerical models” section, the prediction accuracy of the bubble behaviors and phase distribution characteristics was mainly affected by the model of interphase momentum transfer, including the model of drag force, lift force, wall lubrication force, turbulence dispersion force, etc. The Ishii Zuber drag force model was adopted for all the cases studied in this paper, for which has been proved to be the appropriate model for bubbly flow in round tubes by several authors. Wall lubrication force is the kind of force that tends to push the bubbles away from the wall, which is relevant to the mean diameter and surface tension between the two phases. In this paper three-wall lubrication force model was tested, that is the Antal model, Tomiyama model, and Frank model. To be convenient for the comparison, the void fractions, at sections of 0°, 90°, and 180° adopted these three models, respectively, are shown in Figure 3.

Figure 3.

Comparison of different wall lubrication force model. (V_l=1 m/s, V_g=0.131 m/s). (a)-(c) Antal model; (d)-(f) Tomiyama model; (g)-(i) Frank model; (j)-(l) experiment data by Usui.

It can be seen that the Antal model predicts the void fraction profile of 0° section quite differently from the others because the peak of void fraction appears on the wall, while Tomiyama model and Frank model predict the peak being a distance away from the wall. This is the evidence that Antal model is not able to sufficiently counterbalance the Tomiyama lift force in the near-wall region. Meanwhile, the void fraction peak of 90° section predicted by Antal model and Tomiyama model was farther away from the inner wall of the bend than the experimental result and the simulation results by the other models, which means the distance of the wall lubrication force being active is overestimated. Also the height of void fraction peak of 90° section is underestimated by Antal model and Tomiyama model. These results show that the Frank model may be a better choice to predict the void fraction distribution in the entire region of the U-bend.

Void distribution characteristics under different flow conditions

The local void distribution profile in the cross sections of the U-bend changes continuously with the proceeding of the fluid, as shown in Figure 4. In the straight tube before the inlet, the bubbles concentrate in the narrow region near the wall but do not attach to the wall, which makes the profile to be symmetric about the axis of the tube and similar to that of upward bubbly flow in straight tubes. However, after entrance to the bend, the buoyancy force and centrifugal force began to be active. For the cases with higher liquid superficial velocity (Figure 4(a) to (c)) the centrifugal force makes the liquid flow along the outer side wall and drives the bubbles to move toward the inner wall. A single peak was formed finally in the 90° section, which is attached to the wall closely. This means that in these cases the centrifugal force is too large to be counterbalanced by the buoyancy force and wall lubrication force. After 90° section, the single peak began to leave the inner wall and separate into two peaks. However, the mechanism of this process is still not very clear and needs to be studied further. For the case with lower liquid superficial velocity, the buoyancy force will be dominant and keeps the bubbles flow away from the inner wall, as shown in Figure 4(d).

Figure 4.

Void distributions in symmetry plane, 90° section, and 180° section under different flow conditions. ((a) V_l=4 m/s, V_g=0.131 m/s; (b) V_l=4 m/s, V_g=0.5 m/s; (c) V_l=4 m/s, V_g=0.5 m/s; and (d) V_l=4 m/s, V_g=1.0 m/s.)

The void fraction distribution moves away from the inner wall with the increase of water superficial velocity and the decrease of air superficial velocity, because of the increase of the centrifugal force. Also the traces of the bubbles leaving the outer side wall to the inner side wall at the inlet of U-bend, and that of the bubbles attaching the inner wall, were both enlarged with the increase of air superficial velocity.

Conclusions

In this paper the two-fluid model and population balance method were introduced into the simulation of adiabatic water/air bubbly two-phase flow. The interphase momentum transfer process is crucial to the bubble behavior and phase distribution. Through the comparison between simulation and experimental results, we found that Antal wall lubrication model is not able to counterbalance the Tomiyama lift force in the near-wall region. The Antal and Tomiyama model overestimated the active distance of wall lubrication force. However, the Frank model can get the results that meet the experimental results reasonably well if we choose the model parameters carefully. We also simulated the phased distributions under different flow conditions and find that under lower liquid superficial velocity, the buoyancy force is dominant which makes the bubbles concentrate on the outer side of the tube wall; while under higher liquid superficial velocity; the centrifugal force is dominant and makes the bubbles concentrate on the inner side of the tube wall.