Modeling of bubble flows in vertical pipes with the N-phase compressible algebraic slip model

Accepted Manuscript

Modeling of bubble flows in vertical pipes with the N-phase compressible Algebraic Slip Model B.K. Agarwal, C. Narayanan PII: DOI: Reference:

S0301-9322(17)30687-0 10.1016/j.ijmultiphaseflow.2018.04.018 IJMF 2797

To appear in:

International Journal of Multiphase Flow

Received date: Revised date: Accepted date:

13 September 2017 6 March 2018 29 April 2018

Please cite this article as: B.K. Agarwal, C. Narayanan, Modeling of bubble flows in vertical pipes with the N-phase compressible Algebraic Slip Model, International Journal of Multiphase Flow (2018), doi: 10.1016/j.ijmultiphaseflow.2018.04.018

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • N-phase compressible algebraic slip model (ASM) for bubbly flows in vertical pipes.

CR IP T

• Gas number density should be solved to account for full effects of compressibility. • Simulation results are able to predict the behavior of different types of bubbly flow.

• Inadequacy of the current wall lubrication force models is clearly shown.

AC

CE

PT

ED

M

AN US

• The ASM model is shown to be as accurate but computationally efficient compared to two-fluid models such as MUSIG model.

1

ACCEPTED MANUSCRIPT

Modeling of bubble flows in vertical pipes with the N-phase compressible Algebraic Slip Model

a

CR IP T

B. K. Agarwala , C. Narayananb,∗

Eidgen¨ ossische Technische Hochschule Z¨ urich, Department of Mechanical and Process Engineering, ETH-Zentrum, 8092 Z¨ urich, Switzerland b ASCOMP AG, Technoparkstrasse 1, 8005 Zurich, Switzerland

AN US

Abstract

AC

CE

PT

ED

M

In the current work, the N-phase compressible algebraic slip model (ASM) c which framework is presented. The framework is implemented in TransAT is a commercial finite-volume CFD solver specialized in the modeling of multiphase flows. It is applied to bubbly flows in vertical pipes where compressibility of the gas phase is important. Effects of compressibility of the gas phase namely, the change in volume of the individual bubbles, the reduction in the liquid volume fraction, and the change in bubble size due to almost constant number density are identified and illustrated with the help of an analytical one-dimensional model. It was shown that the volume effect was the strongest, although the effect of the change in the bubble radius is also important though to a lesser extent. c does not solve for the gas The framework implemented in TransAT number density and assumes a fixed radius for the bubbles. In the first step the implementation is verified with the help of the one-dimensional model. In the next step, the implementation is tested against the selected test cases from the TOPFLOW database for different types of bubbly flow namely the wall-peaked, the transition and the central-peaked. The simulation results are able to predict the overall behavior of the different types of bubbly flow. The radial gas velocity profiles for all the test cases showed a good match with the experimental data. Unrealistic accumulation of small bubbles close to the walls is observed for the wall-peaked test cases. The results showed that the unrealistic behaviour could be removed tuning the wall lubrication force model. The results highlight the inadequacy of the current wall lubrication force models and show that the wall lubrication force should depend ∗

Corresponding author Email address: [email protected] (C. Narayanan)

Preprint submitted to International Journal of Multiphase Flow

April 30, 2018

ACCEPTED MANUSCRIPT

on the velocity of the continuous phase probably in an analogous way to the lift force.

CR IP T

Keywords: compressible multiphase flow, bubbly flow, mixture model, number density, wall lubrication force, TOPFLOW 1. Introduction

AC

CE

PT

ED

M

AN US

Multiphase flows can be found in various industries like nuclear, chemical, oil and gas etc. Over the last few decades, CFD has been used widely for predicting the behavior of such flows. Accurate prediction helps in improving the efficiency of the industrial processes and also contributes to the safety aspect especially in the nuclear industry. Multiphase flows can be broadly classified into three types, separated flows (a distinct interface exists between two-phases e.g., annular flows, stratified flows), dispersed flows (bubbly flows, solid particles in liquid, droplets or solid particles in gas) and transitional flows. Bubbly flow which forms a part of the dispersed flows occurs when the gas volume fractions ranges from low to intermediate level. In certain applications concerning bubbly flows, it is found that compressibility of gas phase is very important. E.g., in bubbly flows in long vertical pipes, the change in hydrostatic head of the liquid phase introduces a change in the density of the gas phase. Due to the conservation of gas mass flux, a corresponding change takes place in its volume flux. This change in the volume flux of the gas phase causes a change in its volume fraction and the diameter of the bubble. Hence, modeling of these kinds of flows ignoring compressibility of the gas phase leads to under prediction in the gas volume fraction. In the literature, bubbly flows have been by mainly dealt with the MUSIG model (Lo, 1998) under the Eulerian framework. The homogeneous MUSIG model takes into account a poly-dispersed bubble distribution size by solving mass conservation equation for multiple bubble size groups. It also solves for mass and momentum conservation equation for the liquid phase while only one momentum equation for the entire group of bubbles is solved. This limits the applicability of this model to applications where the difference between the velocities of different bubble groups is not important. It is known that coefficient of lift force acting on bubbles changes its sign based on bubble size. As such, this model cannot predict radial de-mixing of bubbles as has been reported in Krepper and Prasser (2000). The inhomogeneous MUSIG model (Shi et al., 2004) addresses this problem by solving 3

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

momentum conservation equations for each bubble size group (Frank et al., 2008). The above method however, is computationally very expensive since the momentum equations have to be solved for each bubble size group. A simpler approach involves solving the governing equations only for the fluid mixture (Manninen and Taivassalo, 1996). The influence of the individual phases are taken into account by means of an algebraic slip model (ASM). This model derives an algebraic expression for the slip velocity between the phases by a balance of the interfacial forces and the pressure gradient. This step is same as the modeling of the interfacial term in case of the MUSIG model. However, because of the model’s basic assumption of instant relaxation of the gas phase velocity to its terminal velocity, virtual mass and other transient effects cannot be accounted for in this framework. In general, the time taken by bubbles to reach the terminal velocity as compared to characteristic flow time scales is very small but in certain regions such as gas injection systems, this is not true and the ASM model is not applicable. In literature, the ASM model in combination with a population balance model have been successfully used for solution of bubbly flows (Swiderski et al. (2016) and Bhole et al. (2008)). Thus, we see that within the range of its applicability, ASM model is a simpler and computationally efficient approach for the solution of bubbly flows. Another advantage of this model is its capability to model a poly-dispersed bubble distribution with individual velocities for multiple bubble size groups with very limited additional computational resources. The current work is an extension of the 4+ -equation compressible multiphase framework presented by Labois and Narayanan (2017) by introducing the algebraic slip model (ASM) as proposed by Manninen and Taivassalo (1996). The framework is discussed with the help of an one-dimensional analytical model. The framework without the number density is implemented c a commercial finite volume CFD solver specialized in the in TransAT , modeling of multiphase flows and applied to the experimental database of TOPFLOW (Beyer et al., 2008) for a range of test cases. In this study the influence of gas-phase compressibility is analyzed in detail and exemplified with an experimental setup where this is important (long pipe such that the pressure drop due to gravity and shear will cause a significant change in the gas density). The different factors that influence slip due to the introduction of gas-phase compressibility - such as the change in the density, change in volume fraction and number density have been quantified for a simple 1-D problem. The importance of solving a number density equation for compressible phases has been brought out. A path to 4

ACCEPTED MANUSCRIPT

CR IP T

improve the wall lubrication force model currently in use has been proposed. The study shows that the ASM model can deliver results of the same or better quality compared to MUSIG model due to the fact that the basic limitation for bubbly flow models lies in the inadequacy of the interphase force models, such as the wall lubrication model. 2. Mathematical Formulation

AN US

equations derived from the conservation laws. It involves solution of a continuity equation for each phase while the momentum, the temperature and the pressure equations are solved for the fluid mixture. Additionally, a number density equation may be solved for phases which are compressible. The slip velocity between the phases is calculated algebraically.

2.1. Governing equations The governing equations for the N-phase compressible multiphase flows are as follows.

M

2.1.1. Phase mass conservation equation The mass conservation equation for a phase k in terms of the volume fraction (αk ) is written as,   ∂ ∂  (αk ρk ) + αk ρk uj + upd ˜td = j jk + u ∂t ∂xj

νt ∂αk ρk Sct ∂xj

ED

∂ ∂xj

!

000

+m ˙k

(1)

AC

CE

PT

where ρk is the density given by an equation of state, uj is the mass-weighted mixture velocity, upd jk is the drift velocity due to pressure gradient (Manninen 000 and Taivassalo, 1996) and m ˙ k is the mass transfer rate per unit volume between the phases. νt = µt/ρ where µt is the turbulent eddy viscosity and ρ the density of the fluid mixture. Sct is a turbulent Schmidt number. The other terms in the equation represent the contribution of turbulence drift based on the turbulent slip velocity as defined in Manninen and Taivassalo (1996). 2.1.2. Phase number density equation For a mono-disperse bubbly flow, the volume fraction of a phase k can be expressed as the product of the phase number density (nk ) and the volume of an individual bubble (vbub ) as, αk = nk vbub 5

(2)

ACCEPTED MANUSCRIPT

A transport equation for nk without phase change, bubble coalescence and breakup effects can be expressed as (Swiderski et al., 2015),

2.1.3. Mixture momentum equation The mixture momentum equation is written as, ∂ ∂p ∂ ∂ (ρui uj ) = − + (ρui ) + ∂t ∂xj ∂xi ∂xj N X

∂ − ∂xj

αk ρk udik udjk

+ M σ + ρgi

!

(3)

2 2 (µm + µt ) (Sij − Skk δij ) 3



(4)

AN US

k=1

!



νt ∂nk Sct ∂xj

CR IP T

 ∂ ∂ ∂   nk uj + upd ˜td = (nk ) + j jk + u ∂t ∂xj ∂xj

where one of the key assumptions is the existence of a mixture viscosity model (rheology) for µm that accounts for the viscous momentum transfer in the fluid mixture. udik is the total drift velocity due to pressure gradient and turbulence, p is the mixture pressure, gi is the gravity vector, and Sij is the symmetric strain-rate tensor. M σ represents the effect of surface tension which is neglected in the current modeling (Labois and Narayanan, 2017).

M

2.1.4. Mixture temperature equation The mixture temperature equation (Labois and Narayanan, 2017) can be written as, !

ED

ρCp

∂T ∂T + uj ∂t ∂xj 000

PT

+Φ + Q +

P X

"

∂p ∂ ∂p ∂T + uj + = (λ + λt ) ∂t ∂xj ∂xj ∂xj

000

m ˙ n Ln

#

(5)

n=1

AC

CE

where Cp is the mixture sensible enthalpy, Φ is the viscous dissipation of kinetic energy and the last term is the latent heat contribution due to P different components undergoing phase change. The above equation can be simplified by writing in terms of material derivative advected by the mixture velocity for T and p as, 1 Dm T = Dt ρCp



Dm p +Ξ+M Dt



(6)

where Ξ represents the entropic part of heat transfer (diffusion, viscous dissipation and volumetric heating) and M represents the net phase change term of Eq. (5). 6

ACCEPTED MANUSCRIPT

" X αk

ρ

k

−

!#

" # X αk ∂ρk + ρ ρ ∂T k

k

Dm p Dt

Ξ+M ρCp

p

!

  ∂uj X 1 ∂  pd + ρk αk u ˜td j + ujk ∂xj ρk ∂xj k

X 1 k

ρk

"

∂ ∂xj

2.2. Algebraic Slip Model

!

#)

AN US

= −ρ

(

ρk



1 ∂ρk ∂ρk + ∂p T ρCp ∂T p

CR IP T

2.1.5. Mixture pressure equation The mixture pressure equation for compressible multiphase flows has been derived by Labois and Narayanan (2017). The final form is given by,

νt ∂ 000 (αk ρk ) + m ˙k σt ∂xj

(7)

ED

M

Algebraic slip model by Manninen and Taivassalo (1996) involves obtaining an algebraic closure relation for the drift velocity (upd jk ) induced by pressure gradient. This is done by equating the different interfacial forces acting on the dispersed phase with the force due to pressure gradient (buoyancy and/or other effects such as rotation, acceleration, etc). The force balance is written as, Fk = αk

ρ − ρk ∇p ρ

(8)

PT

where Fk is the net force acting which includes forces such as drag, lift, etc.

AC

CE

2.2.1. Drag force The relative motion of bubbles with respect to the continuous phase leads to a drag force acting on the bubbles. The drag force expression from Ishii and Mishima (1984) is given by, αk ρp s s 3 Fkd = − CD ukp ukp 8 Rk

(9)

where ρp is the density of the continuous phase, Rk is the radius of the bubble and CD is the drag coefficient.

7

ACCEPTED MANUSCRIPT





ρp uskp 2Rk

and friction µp factor F, a general expression for the drag coefficient can be written as, Introducing the bubble Reynolds number Reb = 24 F Reb

(10)

CR IP T

CD =

AN US

In the current work, the closure for F has been used from the model derived by Tomiyama et al. (1995). This model takes into account the increase in drag force due to the deformation of bubbles. This deformation is charac|ρp − ρk | |g| (2Rk )2 terized using the E¨ otv¨ os number, Eo = , where g is the σ gravitational acceleration and σ the surface tension. The friction factor is then given by, 

F = max 1 + 0.15Re0.687 , b

1 Eo Reb 9 Eo + 4



(11)

M

2.2.2. Lift force Bubbles moving in a shear flow experience a lift force perpendicular to its motion. The expression of the lift force is given by Drew and Passman (1998). Fkl = −CL αk ρp uskp × (∇ × up ) (12)

CE

PT

ED

where CL is the lift coefficient. Tomiyama et al. (1995) studied the trajectory of single bubbles rising in a simple shear flow of a glycerol water solution and found that the lift coefficient changes its sign for bubbles above a critical diameter. The lift coefficient is positive for small bubbles without any substantial deformation. As such the small bubbles are pushed in the direction of decreasing velocity of the continuous phase. Bubbles above the critical diameter have a negative lift coefficient and hence they are pushed in the direction opposite to that of small bubbles. The correlation for the lift coefficient derived by Tomiyama is given by,

AC

CL =

   min (0.288tanh (0.121 Reb ) , f (Eod )) Eod < 4

4 ≤ Eod < 10 Eod ≥ 10

f (Eo )

d   −0.27

(13)

where f (Eod ) = 0.00105 Eo3d −0.0159 Eo2d −0.0204 Eod +0.474, and Eod , the modified E¨ otv¨ os number defined based on the longest axis dH of a deformed bubble is given by, Eod =

|ρp − ρk | |g| (dH )2 σ 8

(14)

ACCEPTED MANUSCRIPT

The longest axis of a deformed bubble can be obtained from the empirical correlation for the aspect ratio proposed by Wellek et al. (1966) as, = 2Rk

q 3

1 + 0.163Eo0.757

(15)

CR IP T

dH

The variation of CL with radius of bubble for air/water flows is shown in Fig. 1. The sign of CL changes at Rk = 2.9 mm. 4 distinct regimes are identified as per the behavior of CL are,

M

AN US

1. Constant positive regime - In this regime, CL has a constant positive value. This regime ranges from 0.2 ≤ Rk ≤ 2.2 mm. 2. Linearly varying positive regime - In this regime, CL reduces linearly with a constant negative slope and attains a value of 0 at the end of this regime. This regime ranges from 2.2 < Rk ≤ 2.9 mm. 3. Linearly varying negative regime - In this regime, CL is negative. The absolute value of CL increases linearly with a constant negative slope. This regime ranges from 2.9 < Rk ≤ 3.6 mm. 4. Constant negative regime - In this regime, CL has a constant negative value. This regime is valid for Rk > 3.6 mm.

ED

A

B

C

PT

A - Constant positive regime B - Linearly varying positive regime C - Linearly varying negative regime D - Constant negative regime

CE

D

AC

Figure 1: Variation of lift coefficient CL with radius of bubble for air/water flows.

c the lift force expression (Eq. 12) is multiplied by a factor In TransAT , to reduce its value close to the walls. This is done to prevent agglomeration of small bubbles close to the walls. Fkl = −CL αk ρp uskp × (∇ × up ) (1 − exp (−yw/2Rk ))

where yw is the distance of the center of the bubble to the wall. 9

(16)

ACCEPTED MANUSCRIPT

Fkw =

2 1 CW L ρp αk uskp nw Rk

CR IP T

2.2.3. Wall lubrication force Bubbles moving close to the wall in a quiescent liquid experiences a force away from the wall due to the hydrodynamic pressure difference acting on the bubbles. The expression for wall lubrication force derived by Antal et al. (1991) is given by, (17)

where nw is the unit vector perpendicular to the wall, pointing into the fluid domain. CW L is the wall coefficient which depends on Rk and yw . The coefficient is positive and hence bubbles are driven away from the wall. Rzehak et al. (2012) suggested an expression for CW L as, 

AN US

CW L



Rk = max 0, −0.005 + 0.05 yw

(18)

With the above defined forces, the force balance given by Eq. (8) can be rewritten as,

M

Fkd + Fkl + Fkw = αk

ρ − ρk ∇p ρ

(19)

3 αk ρp s s CD ukp ukp − CL αk ρp uskp × (∇ × up ) (1 − exp (−yw/2Rk )) 8 Rk 2 1 CW L ρp αk uskp nw = αk (ρ − ρk ) g (20) Rk

CE

−

PT

ED

∇p is calculated from the mixture momentum equation. In the current implementation, it is approximated by the product of the mixture density and a component of gravity vector ρg. This simplification is justified as the frictional and acceleration pressure drop in cases under consideration are less than 2 %. Plugging the expressions for each force in Eq. (19),

+

AC

Eq. (20) is solved using Newton-Raphson method to obtain the slip velocity. Knowing the slip velocity, the drift velocity due to pressure gradient upd jk can be computed from the following expression, upd jk =

N X

Yp usjkp

p=1

where Yp is the mass fraction of the phase p. 10

(21)

ACCEPTED MANUSCRIPT

2.3. Equation of state Density of each phase is obtained from an equation of state defined for that phase. In the current work, the ideal gas equation given by, p RT

(22)

CR IP T

ρ=

is used for the gas phase, where R is the specific gas constant. For the liquid phase, the stiffened gas equation of state is used. It reads, ρ=

γ(p + p∞ ) Cp (γ − 1)T

(23)

AN US

where Cp and γ are the specific heat capacity at constant pressure and specific heat ratio for the liquid. p∞ is the cohesion pressure of the liquid. 2.4. Compressibility effects The differential of Eq. (2) for the gas phase can be expressed as, dαg = ng dvbub + vbub dng

(24)

PT

ED

M

It is seen from the above equation that the change in gas volume fraction is composed of two effects - 1) change in volume of the individual bubbles and 2) change in the gas number density. Under incompressible conditions, the gas volume fraction changes only due to the second effect. The compressibility of the gas phase changes the volume of the individual bubbles. This is the direct effect of compressibility of the gas phase on its volume fraction. The change in volume of the bubbles also changes the slip velocity which brings about a change in the gas number density. Thus, we see that for a compressible gas phase, both the effects are coupled.

CE

3. Numerical Discretization

AC

c software developed by ASThe model is implemented in TransAT COMP. TransAT is a finite-volume CFD solver specialized in the modeling of multiphase flows. The convective terms in the momentum and phasic mass conservation equations are discretized using the HLPA scheme (Zhu, 1991) whereas in the pressure equation the first-order upwind is employed. The diffusion terms are discretized by second-order central differencing in the momentum and temperature equations. The pressure equation is solved using the SIMPLEC pressure correction method of van Doormaal and Raithby (1984). The simulations presented here are run in steady state and a convergence tolerance of 1e-5 is maintained. 11

ACCEPTED MANUSCRIPT

4. One-dimensional analytical model

4.1. Description of the model

AN US

CR IP T

c The N-phase compressible ASM framework implemented in TransAT does not solve the gas number density equation; it assumes a fixed radius for the bubble. Thus, it ignores the change in drift velocity of the gas phase due to the change in radius of the bubble. The effect of compressibility is reflected only as a net change in the volume fraction of the gas phase. Here, we consider a one-dimensional frictionless pipe in order to understand the full compressible effects. However, in the first step, the gas number c density equation is not solved similar to the implementation in TransAT . This is done to verify the implementation of the N-phase compressible ASM c The focus in this step is to check the correctness framework in TransAT . of the implementation.

ED

M

A mixture of water and air is injected into a 3 m long pipe. The pressure at the inlet of the pipe is 101300 Pa. The volume fractions of the phases at the inlet are: αl = 99 % and αg = 1 %. The temperature of the fluid mixture is set to 293 K. The problem is steady and adiabatic. At the inlet of the pipe, the slip velocity is set to be zero and velocity of the fluid mixture is set to 0.02 m/s. The properties of water and air used in the current analysis are shown in Table 1. The radius of the bubble is set to 60 µm at the inlet boundary. Table 1: The properties of water and air.

Property

PT

Water

AC

CE

Molecular Weight [g/mol] Specific Heat Ratio [−] Cohesion Pressure [Pa] Viscosity [Pa.s] Specific Heat [J/kg.K] Thermal Conductivity [W/m.K]

Value Air

18.015 16.64 1.35e8 9.772e-04 491.01 0.60475

28.97 1.4 1.511e-05 1005 0.0257

For the above described one-dimensional pipe, starting from the governing equations described in Section 2, analytical expressions can be obtained for all the relevant quantities like the mixture pressure, the volume fractions and the drift velocities of the phases. The detailed derivation is presented in Appendix A. 12

ACCEPTED MANUSCRIPT

4.2. Simulation setup

g = 9.81 m/s2 3m =1% = 99 % = 0.02 m/s = 293 K

Y

x

PL = 72184 Pa

AN US

αg αl u T

CR IP T

c A schematic The one-dimensional pipe is simulated using TransAT . of the setup with the inlet and the outlet boundary conditions is shown in Fig. 2. Since the outlet is not available, the outlet pressure computed from the analytical solution is used for the simulation. Velocity is specified at the inlet and at the outlet, the pressure is set.

Figure 2: Simulation setup of the one-dimensional analytical model.

4.3. Results and discussion

AC

CE

PT

ED

M

c output and the anaFigure 3 shows a good match between TransAT lytical results. The axial pressure variation in Fig. 3(a) computed from the c analytical model agrees well with the numerical solution from TransAT which solves the full governing equations described in Section 2. This justifies the assumption made in the analytical solution to neglect the contribution of the viscous term and the drift velocity term in the mixture momentum equation. The mixture density changes by a small amount due to a small change in the gas volume fraction as seen in Fig. 3(b). Thus, the axial pressure distribution is largely driven by the change in hydrostatic head of the liquid. A sharp drop in the gas volume fraction close to the inlet boundary is seen in Fig. 3(b). At the inlet, the slip velocity is set to zero. However, the slip velocity due to the ASM model is fully applied at the downstream cell face of the first computational cell. This difference in the slip velocities at the inlet and the downstream face of the first computation cell is the reason for the observed drop in the gas volume fraction close to the inlet. After the first cell, a gradual increase in the gas volume fraction is seen. This is due to the decrease in the density of the gas phase. The density of the gas phase decreases due to the decrease in hydrostatic pressure of the liquid. Since the slip velocity is nearly constant, the conservation of the gas phase mass flux ensures a corresponding increase in the gas volume fraction. 13

ACCEPTED MANUSCRIPT

0.105

1.05

(a)

Analytical - Fixed radius TransAT - Fixed radius

(b)

1.00 0.95

0.090

0.90

αg [%]

0.095

0.085 0.080

Analytical - Fixed radius TransAT - Fixed radius

0.85 0.80

0.075

0.75

0.070 0.0

0.70 0.0

0.5

1.0

1.5 x [m] 7.96

2.0

2.5

(c)

udg × 1000 [m/s]

7.94

7.93

M

7.92

7.90 0.0

0.5

1.0

1.5 x [m]

2.0

2.5

3.0

Analytical - Fixed radius TransAT - Fixed radius

7.95

7.91

3.0

AN US

Pressure[MPa]

0.100

CR IP T

Fig. 3(c) shows that the slip velocity of the gas phase decreases axially due to the dependence of the slip velocity on the liquid volume fraction and the assumption of a fixed radius for the bubble.

0.5

1.0

1.5 x [m]

2.0

2.5

3.0

ED

c output (a) pressure, (b) gas volume fraction and (c) Figure 3: Comparison of TransAT gas drift velocity with results from the analytical model for fixed radius of the bubble.

AC

CE

PT

In the second step, for the analytical model, the gas number density is also solved in addition to the gas volume fraction. This enables the consideration of the full compressibility effects of the gas phase. As seen from Fig. 4(a), the slip velocity increases axially contrary to the previous result. Thus, the effect of the increase in radius of the bubbles on the slip velocity is much more dominant than the effect of the decrease in the liquid volume fraction. The new gas volume fraction computed is plotted along with the results from the fixed radius assumption in Fig. 4(b). The individual contributions of the change in volume of the individual bubbles and the change in the gas number density as discussed in Section 2.4 are also plotted. The lower value of the gas volume fraction is due to the increase in the slip velocity. c and the analytical model with the assumption of At the outlet, TransAT the fixed radius over predicts the gas volume fraction by ∼ 7% as compared 14

ACCEPTED MANUSCRIPT

0.010

CR IP T

to the analytical model with the gas number density. Thus, for an accurate prediction, the number density should be solved in addition to the volume fraction for all the compressible phases. (a)

0.006 0.004

0.000 0.0

1.05 1.00 0.95

1.0

1.5 x [m]

2.0

2.5

3.0

2.5

3.0

Analytical - Fixed radius TransAT - Fixed radius Analytical - Variable radius Vol. change contribution ng change contribution

0.85

ED

αg [%]

0.90

(b)

0.5

AN US

Analytical - Fixed radius TransAT - Fixed radius Analytical - Variable radius

0.002

M

udg [m/s]

0.008

0.80

0.75

PT

0.70

0.65 0.0

0.5

1.0

1.5 x [m]

2.0

AC

CE

Figure 4: Comparison of results (a) drift velocity of gas (b) gas volume fraction from the analytical model for fixed and variable radius. For the variable radius, contribution of the change in volume of the bubble and the change in the gas number density are also plotted.

5. Application to the TOPFLOW experiment c The N-phase compressible ASM framework implemented in TransAT is applied to the TOPFLOW (Beyer et al., 2008) results which is one of the most comprehensive experimental databases concerning two-phase flow 15

ACCEPTED MANUSCRIPT

that are available today. A brief description of the experimental setup, the selection procedure of the test cases, the simulation setup and the results are presented next.

CR IP T

5.1. Description of the experimental setup

M

AN US

A vertical pipe of inner diameter of Dpipe = 195.3 mm and an approximate height of L = 9 m was used for the investigation of air-water flows. A constant temperature of 300 C ± 1 C was maintained. A two-level low temperature wire mesh sensor was mounted at the top of the vertical pipe whose position cannot be altered. In order to be able to study evolution of gas along the height of the pipe, injected at a particular location, a variable gas injection module capable of injecting at multiple locations was attached to the vertical pipe. A constant gas injection pressure of 0.25 MPa was maintained at each injection location. Thus, with these arrangements, the measurements carried out at the sensor location for different injection locations is equivalent to having a fixed injection location and measurements being carried out at different locations along the pipe. The variable gas injection module allows injection of gas from six different locations distributed over the total height of the pipe. Gas is injected through the wall orifices. The distances of the gas injection locations from the measurement plane are typically expressed in terms of L/Dpipe ratio and are listed in Table 2.

ED

Table 2: Gas injection length from the measurement plane expression in mm and L/Dpipe ratio.

A

I

L

O

R

Length, L (mm) L/Dpipe ratio

221 1.1

1552 7.9

2595 13.3

4531 23.2

7802 39.9

PT

Level

AC

CE

The database provides time averaged radial gas volume fraction profiles, radial distribution of diameter of bubbles and radial gas velocity profiles. It has been pointed by Rzehak and Krepper (2015) that the mass fluxes computed from the experimental gas volume fraction and gas velocity profiles are higher than the inlet flow parameters. In the bubbly flow regime, they found the mass fluxes to be approximately higher by a factor of 1.2 which they used to divide the experimental gas volume fractions in their work. Throughout this work, the experimental radial gas volume fraction profiles used for comparison are also divided by the same factor.

16

ACCEPTED MANUSCRIPT

5.2. Selection of the test cases

PT

ED

M

AN US

CR IP T

The current work focuses on the application of the ASM model for bubbly flows. It assumes that bubble coalescence and breakup effects are negligible which is true in the case of the TOPFLOW experiments considered in this study (Beyer et al., 2008). It has been pointed out in the report by Beyer et al. (2008) that based on the radial gas volume fraction distributions at the maximum length of the pipe, bubbly flows can be broadly categorized into three different types. If the radial distribution shows a maximum close to the walls, it is termed as wall-peaked and that peaking at the center of the pipe is termed as central-peaked. The profiles in between these two types are termed as transition. The three different types of bubbly flow for a fixed jl = 1.017 m/s and varying jg as observed in the TOPFLOW tests are shown in Fig. 5.

CE

Figure 5: Radial gas volume fraction distribution at level R (L/Dpipe = 39.9) with marked bubbly flow types. jl is kept constant at 1.017 m/s while jg is varied from 0.0096 to 0.0368 m/s (Beyer et al., 2008).

AC

Two test cases from each type are selected in the current study such that they have different flow rates of water, except for the wall-peaked bubbly flow as no wall-peaked bubbly flow was observed with lower flow rate for water for the given flow rates of air. An overview of the selected test cases is presented in Table 3. The superficial velocities (jl , jg ) and Reynolds number are tabulated for each test case.

17

ACCEPTED MANUSCRIPT

Table 3: Superficial velocities (jl , jg ) and Reynolds number for the selected test cases.

Case jl (m/s) jg (m/s) Re

008 1.017 0.0025 202803

Transition

041 1.017 0.0096 202803

006 0.405 0.0025 80762

052 1.017 0.0151 202803

Central-peaked 039 0.405 0.0096 80762

074 1.017 0.0368 202803

CR IP T

Wall-peaked

AN US

5.3. Simulation setup c The three-dimensional The test cases are simulated using TransAT . vertical pipe is modeled as an axisymmetric case (as by Rzehak and Krepper (2015)). A schematic of the setup is shown in the Fig. 6. Wall

ED

M

Inflow

R = 97.65 mm

Level R

Outflow

g = 9.81 m/s2

Level A

203 mm

7581 mm

PT

Axis of symmetry Figure 6: A schematic of the setup used for the simulation of TOPFLOW cases.

AC

CE

At the inflow boundary, the phasic velocities and the phasic volume fractions are specified. In order to simulate injection through the wall orifices, injection was carried out through a small annular ring close to the walls. This methodology however, introduced a dependence of results on the inlet boundary conditions. Hence, it was discarded and the experimental radial volume fraction profiles at level A are used as the inlet boundary conditions. Constant phasic velocities (uk ) are specified at the inlet boundary so as to match the inlet phasic volume fluxes as obtained from Eq. (25). uk =

π (D2pipe/4) jk 2π

R Dpipe/2 0

αk rdr

(25)

The outlet boundary is located at around 200 mm from level R so as to 18

ACCEPTED MANUSCRIPT

pout = pinj − ∆p

CR IP T

nullify any effect of the outlet boundary condition on the results. At the outlet boundary, pressure (pout ) is specified. This pressure is obtained from Eq. (26). pinj is the gas injection pressure which is equal to 0.25 MPa and ∆p, the pressure drop in the vertical pipe between level A and the outlet boundary is obtained from the experimental results. (26)

A mixture of air and water is injected with a constant temperature of C at the inlet boundary. The properties of water and air used for the simulation are same as those defined in Table 1. The pipe wall is modeled as adiabatic and no-slip boundary conditions are applied for the fluid mixture at the wall. Turbulence is modeled using the RANS approach. The high Re k −  two-equation model with standard wall functions is employed. The ASM model requires the specification of a bubble radius for each phase for the computation of the slip velocity. In Eq. (20), the terms on the left hand side are the different forces acting on the bubble which depend on the shape of the bubble. The term on the right hand side is the buoyancy force acting due to gravity which is a volumetric force. Therefore, Sauter mean diameter is chosen as a representative diameter as it represents a ratio of the bubble volume to the bubble surface area. Sauter mean diameter for a phase k (Dk32 ) is obtained from the experimental bubble size distribution at level A using Eq. (27).

PT

ED

M

AN US

300

R∞ 3 D nk (D)dD R Dk32 = 0∞ 2 0

D nk (D)dD

(27)

AC

CE

where nk (D) represents the number density of the phase k for a polydisperse bubbly flow. Two approaches are used to simulate the test cases. In the first approach, the polydisperse bubble size distribution is approximated by a single representative radius (Rk32 ). This reduces the N-phase ASM framework to 2-phase ASM framework. The advantage of this approach is that it is computationally faster. In the second approach, the distribution is captured by multiple bubble size bins. In the current work, four bins are used with each bin corresponding to the regimes shown in Fig. 1. The advantage of this approach is that each bin has its own velocity which will respond to the different forces acting on them, particularly the lift force more accurately. Rk32 for the single bin and the multiple bins approaches are shown in Table 4.

19

ACCEPTED MANUSCRIPT

Table 4: Sauter mean radius for the single bin and the multiple bins approach.

Single bin

008 041 006 052 039 074

1.92 2.40 2.77 2.84 3.13 4.70

Bin 3 3.71 3.77 3.91 3.90 3.94 4.79

5.4. Mesh dependency study

AN US

CR IP T

Case

Rk32 (mm) Multiple bins Bin 0 Bin 1 Bin 2 1.52 2.34 3.09 1.67 2.42 3.11 1.71 2.50 3.14 1.78 2.51 3.14 1.66 2.55 3.18 1.36 2.51 3.30

AC

CE

PT

ED

M

Mesh dependency study involves increasing the number of the elements of the discretized domain with an aim to obtain grid sizes for which the truncation errors arising from the numerical approximation of the governing differential equations are negligible. With the use of the wall functions for turbulence models, caution has to be exercised so that the reduction in grid sizes is carried out within the validity of the wall function model. The c is valid if y + at standard wall function model implemented in TransAT + the cell adjacent to the wall falls in the range 11.6 ≤ y ≤ 200 − 500. Two different values of jl for the selected test cases warrants an independent study for each jl as y + at the first cell next to the wall depends on the Reynolds number. Starting from a coarser mesh, the total number of elements is increased by a factor of two. The elements in the radial direction are increased by a factor of 22/3 while the elements in the axial direction are increased by a factor of 21/3 as carried out by Frank et al. (2008). In all the cases y + > 11.6 is maintained by adjusting the thickness of the first element next to the wall. The details of the mesh dependency study with the mesh numbers, the number of elements in the radial and the axial direction, the total number of elements, the thickness of the first element and range of y + value for jl = 0.405 m/s and jl = 1.017 m/s are listed in Table 5 and Table 6 respectively. A comparison of the radial gas volume fraction profiles for different mesh numbers with the experimental results at level R for the test cases 039 and 074 are shown in Fig. 7. As expected, the results converge with the increase in number of elements. In Fig. 7(a), a small difference in the results is observed between the mesh number 1 and the mesh number 2. No visible differences are observed between the mesh number 2 and the mesh number 3. 20

ACCEPTED MANUSCRIPT

Table 5: Parameters of the mesh dependency study for jl = 0.405 m/s.

Number of elements in axial direction

Total elements

First element thickness (mm)

1 2 3

28 45 65

778 973 1216

21784 43785 79040

2.8 1.5 1.5

min y +

max y +

CR IP T

Number of elements in radial direction

24 14 14

AN US

Mesh number

50 26 26

Table 6: Parameters of the mesh dependency study for jl = 1.017 m/s.

Number of elements in radial direction

Number of elements in axial direction

1 2 3 4

28 45 72 115

778 973 1216 1520

Total elements

ED

M

Mesh number

21784 43785 87552 174800

First element thickness (mm)

min y +

max y +

1.20 0.75 0.76 0.75

32 20 17 17

51 33 29 29

CE

PT

Hence, mesh independence for jl = 0.405 m/s is obtained with mesh number 2. Similarly, as seen in Fig. 7(b), mesh independence for jl = 1.017 m/s is obtained with mesh number 3. Mesh number 2 and mesh number 3 for jl = 0.405 m/s and jl = 1.017 m/s respectively are used for all further simulations.

AC

5.5. Inlet turbulence parameters The two equation k −  turbulence model used in the current study requires turbulence intensity (TI) and eddy viscosity ratio (EVR) to be specified at the inlet boundary. TI is a measure of turbulent velocity scale. The relation between TI, EVR and the integral turbulent scale lt is given by the equation below. lt = f ×

EVR TI

21

(28)

ACCEPTED MANUSCRIPT

7

(a)

2.5

5

αg [%]

αg [%]

2.0 1.5 1.0 0.5 0.0 0.0

(b)

6

Mesh 1 Mesh 2 Mesh 3 Expt.(Corected)

0.2

0.4 0.6 r/R [−]

4 3 2 1

0.8

1.0

0 0.0

CR IP T

3.0

Mesh 1 Mesh 2 Mesh 3 Mesh 4 Expt.(Corrected)

0.2

0.4 0.6 r/R [−]

0.8

1.0

AN US

Figure 7: Comparison of the radial gas volume fraction profiles for different mesh numbers with experimental (corrected) results at level R for the test cases (a) 039 (jl = 0.405 m/s) and (a) 074 (jl = 1.017 m/s).

ED

M

where f is a factor dependent on the fluid and the flow parameters. In the absence of experimental data for the turbulence quantities at the inlet boundary, a preliminary test was carried with TI = 3.4 % and EVR = 320. The results showed high diffusion near the center of the pipe. Hence, in order to arrive at reasonable inlet turbulence parameters, a parametric study is conducted for different TI and EVR as shown in Table 7. Diffusion could be adjusted by reducing both turbulent velocity (TI) and length scale (lt ). This could be achieved by reducing EVR by a factor larger than the factor used to reduce TI. The factors are chosen such that lt is reduced by a factor of 4.6 in each subsequent test. A similar study of inlet turbulence parameters has been reported in Zschaeck et al. (2014).

Parametric study

TI(%)

EVR

lt

1 2 3 4

10 4.64 2.15 1

1000 100 10 1

100 f 21.6f 4.6 f 1f

AC

CE

PT

Table 7: Parameters of the inlet turbulence study.

The parametric study is carried out for all the selected test cases. The parametric study 3 (TI = 2.15 % and EVR = 10) is chosen for further simulations as it showed a more closer agreement to the experimental results for all the test cases. It was also observed that with the development of the flow, the effect of the inlet turbulence parameters reduces drastically with almost negligible effect at level R. This can be seen from Fig. 8 where a 22

ACCEPTED MANUSCRIPT

(a)

2.5

αg [%]

2.0

2.5 EVR = 1, TI = 1 % EVR = 10, TI = 2.15 % EVR = 100, TI = 4.64 % EVR = 1000, TI = 10 % Expt. (Corrected)

1.5 1.0

1.5 1.0 0.5

0.5 0.2

0.4 0.6 r/R [−]

0.8

1.0

0.0 0.0

EVR = 1, TI = 1 % EVR = 10, TI = 2.15 % EVR = 100, TI = 4.64 % EVR = 1000, TI = 10 % Expt. (Corrected)

0.2

0.4 0.6 r/R [−]

0.8

1.0

AN US

0.0 0.0

(b)

2.0

αg [%]

3.0

CR IP T

comparison of the radial gas volume fraction profiles for different TI and EVR with the experimental results at levels I and O for one of the test case 039 is shown.

Figure 8: Comparison of the radial gas volume fraction profiles for different inlet turbulence parameters with experimental (corrected) results for the test case 039 at (a) level I and (b) level O.

5.6. Simulation results

ED

M

The N-phase compressible ASM framework is employed to simulate the selected test cases using the single bin and the multiple bins approach. The results along with the main observations and discussion is presented next.

AC

CE

PT

5.6.1. Types of bubbly flow Both the approaches are able to predict types of bubbly flow as described in Section 5.2 for each test case (008 and 041 - wall-peaked, 006 and 052 transition and 039 and 074 - central-peaked) accurately with regard to the experimental data. This is depicted in Fig. 9 where the radial gas volume fraction and velocity profiles at the maximum length of the pipe for one test case of each type are shown. The high peak in the gas volume fraction profile at the wall for the wallpeaked case 041 (Fig 9(a)) is unrealistic. The reason for such behavior is discussed in Section 5.6.3. The prediction is reasonably good for the transition and the central-peaked cases. However, during the evolution phase (level I, L and O, not shown here), there is room for improvement for numerical agreement with the experimental data for all the test cases. The radial velocity profiles show a good match to the experimental data at all levels. This points to the applicability of the turbulence model used for the current study. 23

ACCEPTED MANUSCRIPT

1.6

(a)

Expt.(Corrected) Single bin Multiple bins (total)

12

1.2

8

ug [m/s]

αg [%]

10

6 4

0.8 0.6

Expt. Single bin Multiple bins (avg.)

0.2

0 0.0 0.7

0.2

0.4 0.6 r/R [−]

0.8

0.0 0.0

1.0

0.9

(c)

(d)

0.7

0.5 ug [m/s]

0.4

0.4 0.6 r/R [−]

AN US

0.6

0.3

0.8

1.0

0.8

1.0

0.8

1.0

0.5

0.4

0.3

0.1

Expt.(Corrected) Single bin Multiple bins (total)

0.0 0.0

0.2

0.4 0.6 r/R [−]

0.8

M

2.0

1.0

0.2

0.4 0.6 r/R [−]

0.8

0.2

0.4 0.6 r/R [−]

(f)

0.8

0.7 0.6 0.5 0.4 0.3

Expt.(Corrected) Single bin Multiple bins (total)

PT

0.0 0.0

ED

1.5

0.5

0.0 0.0

1.0

0.9

(e)

2.5

Expt. Single bin Multiple bins (avg.)

0.2

0.1

ug [m/s]

0.2

3.0

0.2

0.8

0.6

αg [%]

1.0

0.4

2

αg [%]

(b)

1.4

CR IP T

14

0.2 0.1 1.0

0.0 0.0

Expt. Single bin Multiple bins (avg.)

0.2

0.4 0.6 r/R [−]

CE

Figure 9: Comparison of the radial gas volume fraction (left hand side) and velocity (right hand side) profiles obtained from simulations with experimental results at level R for the test case (a) , (b) 041 (c) , (d) 006 and (e), (f) 039.

AC

5.6.2. Single bin vs Multiple bin In Fig. 9, the sum of the radial gas volume fraction profiles and the massweighted average of the radial velocity profiles of the individual bins in the multiple bins approach are also plotted. Ignoring the unrealistic result of the wall-peaked case 041, the results computed from the single bin approach show only small differences compared to that computed from the multiple bins approach. This suggests that the single bin approach using Sauter mean radius as the representative radius of the bubble is applicable for the 24

ACCEPTED MANUSCRIPT

selected test cases. However, in a more general sense, noting the behavior of the lift force, the multiple bins approach might prove necessary.

AC

CE

PT

ED

M

AN US

CR IP T

5.6.3. Discussion on wall lubrication force The unrealistic accumulation of the bubbles close to the walls for the wall-peaked case 041 (Fig 9(a)) is because of the smaller magnitude of the wall lubrication force as compared to the lift force. This can be solved either by using a lift force model with a lower coefficient or a wall lubrication force model with a higher coefficient. Frank wall lubrication force model (Frank et al., 2008) gives a higher value of coefficient as compared to that used in the current work and so it is also used to simulate the case. Fig. 10 shows that the Frank wall lubrication force model is able to counteract the lift force. Results computed using the inhomogeneous MUSIG model in CFX extracted from Rzehak and Krepper (2015) are also plotted. It can be seen that the peak value of the gas volume fraction profiles computed from the simulations (TransAT and CFX) do not diffuse into the domain with increasing height contrary to the experimental data. This causes differences between the simulation and the experimental results. The multiple bins approach for the transition and the central-peaked cases offers an interesting insight. As seen in Fig. 11, the wall lubrication force is able to counteract the lift force in bin 0 and 1 only for the cases with lower jl which suggests that the lift force is higher for cases with higher jl . A close look at Eq. 16 reveals that the lift force depends on the velocity gradients of the continuous phase. This explains the higher magnitude of the lift force with higher jl . Thus, the wall lubrication force whose form does not depend on the magnitude of the velocity gradients of the continuous phase will be able to counteract the lift force only for a certain range of bulk velocities of the continuous phase and not beyond that. The form of the wall lubrication force used in the current work was derived by Antal et al. (1991) for laminar flows where the velocity of the continuous phase close to the walls was assumed to be small as compared to the slip velocity between the phases and was neglected. Nguyen et al. (2013) showed that the wall lubrication force for turbulent bubbly flows depends on the velocity of the continuous phase. They proposed new correlations for the wall lubrication force coefficients with a dependence on the inlet superficial velocity of the continuous phase. These correlations are also dependent on the pipe diameter. The correlations could not be used for the current work as they depend on global parameters and not on local quantities rendering them inappropriate for a CFD code. Indeed since the time of submission of this article, a publication (Lubchenko 25

ACCEPTED MANUSCRIPT

12

(a)

Antal wall lubrication force Frank wall lubrication force Expt.(Corrected) CFX (Inhom. MUSIG)

10

αg [%]

6

4

4

2

2

12

0.2

0.4 0.6 r/R [−]

0.8

0 0.0

1.0

12

(c)

Antal wall lubrication force Frank wall lubrication force Expt.(Corrected) CFX (Inhom. MUSIG)

10

6 4

0.2 (d)

0.4 0.6 r/R [−]

0.8

1.0

Antal wall lubrication force Frank wall lubrication force Expt.(Corrected) CFX (Inhom. MUSIG)

10 8

αg [%]

8 αg [%]

8

6

0 0.0

Antal wall lubrication force Frank wall lubrication force Expt.(Corrected) CFX (Inhom. MUSIG)

10

AN US

αg [%]

8

(b)

CR IP T

12

6 4

2

2

0.2

0.4 0.6 r/R [−]

0.8

1.0

M

0 0.0

0 0.0

0.2

0.4 0.6 r/R [−]

0.8

1.0

ED

Figure 10: Comparison of the radial gas volume fraction profiles for two wall lubrication force models with experimental (corrected) and CFX results for the test case 041 (wallpeaked) at (a) level I (b) level L (c) level O and (d) level R.

PT

et al., 2018) has appeared that seems to address this issue considering similar arguments to motivate the development of an alternative wall-interaction model.

CE

6. Summary and conclusions

AC

In the current work, an N-phase compressible ASM model framework is presented. It is applied to bubbly flows in vertical pipes where compressibility of the gas phase is important. The two compressibility effects namely, the change in volume of the individual bubbles and the change in gas number density were exemplified with the help of an analytical one-dimensional model. c does not solve for the gas The framework implemented in TransAT number density and assumes a fixed radius for the bubbles. In the first step the implementation was verified with the help of the one-dimensional model.

26

ACCEPTED MANUSCRIPT

αg [%]

0.8 0.6

6 Expt.(Corrected) Bin 0 Bin 1 Bin 2 Bin 3 Multiple bins (total)

5 4

0.4

1

2.5 2.0 1.5

0.4 0.6 r/R [−]

0.8

12 Expt.(Corrected) Bin 0 Bin 1 Bin 2 Bin 3 Multiple bins (total)

10 8

1.0

0.2 (d)

0.4 0.6 r/R [−]

6

0.8

1.0

0.8

1.0

Expt.(Corrected) Bin 0 Bin 1 Bin 2 Bin 3 Multiple bins (total)

4

0.5 0.0 0.0

0 0.0

1.0

αg [%]

(c)

0.2

AN US

0.0 0.0

αg [%]

3

Expt.(Corrected) Bin 0 Bin 1 Bin 2 Bin 3 Multiple bins (total)

2

0.2

3.0

(b)

CR IP T

(a)

αg [%]

1.0

2

0.2

0.4 0.6 r/R [−]

0.8

1.0

0 0.0

0.2

0.4 0.6 r/R [−]

M

Figure 11: Comparison of the radial gas volume fraction profiles at level I obtained from simulations for jl = 0.405 m/s (left hand side) with jl = 1.017 m/s (right hand side) for test case (a) 006 (b) 052 (c) 039 and (d) 074.

AC

CE

PT

ED

A close match was observed between the results indicating the correctness of the implementation. In the next step, the model was tested against experimental results from the TOPFLOW database. Two test cases from each type of bubbly flow namely the wall-peaked, the transition and the central-peaked were simulated using one representative radius in the single bin approach and four multiple bubble bins in the multiple bins approach. The range of the multiple bins were chosen such that they fall in the four distinct regimes of the lift force coefficient CL as shown in Fig. 1. Small bubbles in the first and second bin are pushed towards the wall while large bubbles in the third and fourth bin are pushed towards the center of the pipe. Both the approaches were able to predict the overall behavior of the different types of bubbly flow. With regard to accuracy in the prediction of the radial gas volume fraction profiles during the evolution phase, there is ample room for improvement. The radial gas velocity profiles however, match well with the experimental profiles. Unphysical accumulation of bubbles close to the walls for the wall-peaked cases were predicted. It was also shown that 27

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN US

CR IP T

such behavior could be removed by using a wall lubrication force model with higher coefficient. The results from the multiple bins for the transition and central-peaked cases showed that the wall lubrication force is not able to counteract the lift force in bin 0 and 1 for cases with higher jl . This is because the magnitude of the lift force is higher for higher jl while the opposing wall lubrication force is independent of jl . The non-dependence of the wall lubrication force on the velocity of the continuous phase is an unresolved problem. Nguyen et al. (2013) have pointed out that in the literature, the coefficients of the wall lubrication force has been tuned to match higher lift force arising from the higher velocity of the continuous phase. This tuning of the wall coefficients does not address the root cause of the problem and hence cannot be an effective solution. The work of Nguyen et al. (2013) is a positive step in this direction. Thus, the uncertainties in the modeling of the interfacial forces especially the ratio of the lift and the wall lubrication force coefficient for small bubbles and the non-dependence of the wall lubrication force model on the velocity of the continuous phase should be addressed in the future work. In the analytical one-dimensional model discussed in Section 4, it was seen that ignoring the change in slip velocity due to the change in volume of the bubbles leads to an over-prediction of the gas volume fraction at the outlet of the pipe by ∼ 7%. The parameters of the one-dimensional model such as the ratio of the mixture velocity and the slip velocity and the percentage change in the gas density between the inlet and the outlet boundary are comparable to that of the TOPFLOW cases. As such, we can conclude that a similar over-prediction in the integral of the gas volume fraction at the outlet can result. However, its impact on the radial distribution of the profiles cannot be predicted with the current information. Thus, for a complete representation of the compressible multiphase flows, the gas number density should be solved in addition to the gas volume fraction. References

AC

Antal, S. P., Lahey, R. T., Flaherty, J. E., 1991. Analysis of phase distribution in fully developed laminar bubbly two-phase flow. Int. J. Multiphase Flow 17 (5), 635–652.

Beyer, M., Lucas, D., Kussin, J., Sch¨ utz, P., 2008. Air-water experiments in a vertical dn200-pipe. Technical Report FZD-505, Institute of safety research, Germany: Forschungszentrum Rossendorf. 28

ACCEPTED MANUSCRIPT

Bhole, M. R., Joshi, J. B., Ramkrishna, D., 2008. Cfd simulation of bubble columns incorporating population balance modeling. Chemical Engineering Science 63, 2267–2282.

CR IP T

Drew, D. A., Passman, S. L., 1998. Theory of multicomponent fluids. Springer-Verlag New York, Inc.

Frank, T., Zwart, P. J., Krepper, E., Prasser, H. M., Lucas, D., 2008. Validation of CFD models for mono- and polydisperse airwater two-phase flows in pipes. Nucl. Eng. and Design 238, 647–659.

AN US

Ishii, M., Mishima, K., 1984. Two-fluid model and hydrodynamic constitutive relations. Nucl. Eng. and Des. 82, 107–126. Krepper, E., Prasser, H. M., 2000. Measurements and CFX simulations of a bubbly flow in a vertical pipe. AMIFESF Workshop, Computing Methods in Two-Phase Flow.

M

Labois, M., Narayanan, C., 2017. Non-conservative pressure-based compressible formulation for multiphase flows with heat and mass transfer. Int. J. Multiphase Flow 96, 24–33. Lo, S., 1998. Modeling of bubble breakup and coalescence with the MUSIG model. AEA Technical Report AEAT-4355.

ED

Lubchenko, N., Magolan, B., Sugrue, R., Baglietto, E., 2018. A more fundamental wall lubrication force from turbulent dispersion regularization for multiphase CFD applications. Int. J. Multiphase Flow 98, 36–44.

PT

Manninen, M., Taivassalo, V., 1996. On the mixture model for multiphase flow. VTT Publications, Espoo, Finland.

CE

Nguyen, V. T., Song, C. H., Bae, B. U., Euh, D. J., 2013. The dependence of wall lubrication force on liquid velocity in turbulent bubbly two-phase flows. J. of Nucl. Sc. and Tech. 50, 781–798.

AC

Rzehak, R., Krepper, E., 2015. Bubbly flows with fixed polydiversity: Validation of a baseline closure model. Nucl. Eng. and Design 287, 108–118. Rzehak, R., Krepper, E., Lifante, C., 2012. Comparative study of wall-force models for the simulation of bubbly flows. Nucl. Eng. and Design 253 (0), 41–49.

29

ACCEPTED MANUSCRIPT

CR IP T

Shi, J. M., Zwart, P. J., Frank, T., Rohde, U., Prasser, H. M., 2004. Development of a multiple velocity and multiple size group model for polydispersed multiphase flows. Annual report, Institute of safety research, Germany: Forschungszentrum Rossendorf. Swiderski, K., Caviezel, D., Labois, M., Lakehal, D., Narayanan, C., 2015. Computational modelling of gas-liquid multiphase flows with DQMOM and the N-phase Algebraic Slip Model. In: Computational Methods in Multiphase Flow. Vol. VIII. WIT Press, Valencia, pp. 299–310.

AN US

Swiderski, K., Narayanan, C., Lakehal, D., 2016. Application of N-phase algebraic slip model and direct quadrature method of moments to the simulation of air-water flow in vertical risers and bubble column reactor. Comp. and Chem. Eng. 90, 151–160.

Tomiyama, A., Sou, A., Zun, I., Kanami, N., Sakaguchi, T., 3–7 April 1995. Effects of eotvos number and dimensionless liquid volumetric flux on lateral motion of a bubble in a laminar duct flow. In: Proceedings of the Second International Conference on Multiphase Flow. Kyoto, Japan, pp. 3–16.

M

van Doormaal, J., Raithby, G. D., 1984. Enhancements of the simple method for predicting incompressible fluid flows. Num. Heat Transfer 7, 147–163.

ED

Wellek, R. M., Agrawal, A. K., Skelland, A. H. P., 1966. Shapes of liquid drops moving in liquid media. AIChE J. 12, 854–862.

PT

Zhu, J., 1991. A low diffusive and oscillation-free convection scheme. Comm. Appl. Num. Meth. 7, 225–232.

AC

CE

Zschaeck, G., Frank, T., Burns, A. D., 2014. CFD modelling and validation of wall condensation in the presence of non-condensable gases. Nucl. Eng. and Design 279, 137–146.

30

ACCEPTED MANUSCRIPT

Appendix A. Derivation of the one-dimensional analytical model

  ∂  = 0 ρl αl u + upd l ∂x

CR IP T

The model has been described in Section 4.1. We consider the model to be one-dimensional in the x-direction. The mass conservation equation for the liquid phase is,

(A.1)

Lift and wall lubrication forces are not applicable for one-dimensional model and hence pressure drift is caused only by the buoyancy force acting on the bubbles which is counteracted by the drag force. Considering ρl to be constant, Eq. (A.1) can be rewritten as,

AN US

 ∂   αl u + upd =0 l ∂x

(A.2)

The mass conservation equation for the gas phase is,   ∂  ρg αg u + upd = 0 g ∂x

(A.3)

M

The mixture density can be written as,

ρ = αl ρl + αg ρg

(A.4)

ED

Using the definition of drift velocity for k phases, it can be shown that X

αk ρk upd k =0

(A.5)

PT

k

Summing equations (A.1) and (A.3) and using equations (A.4) and (A.5), the mass conservation equation for the fluid mixture is obtained as,

CE

∂ (ρu) = 0 ∂x

(A.6)

The number density equation for the gas phase can be written as,

AC

 ∂   ng u + upd = 0 g ∂x

(A.7)

Integration of equations (A.2), (A.3), (A.6) and (A.7) between two sections 1 and 2 gives, 



αl u + upd l



2





= αl u + upd l 31



1

= constant = kl

(A.8)

ACCEPTED MANUSCRIPT





ρg αg u + upd g



2





= ρg αg u + upd g



1

= constant = kg

and 



ng u + upd g



2





= ng u + upd g



1

(A.10)

CR IP T

(ρu)2 = (ρu)1 = constant = k

(A.9)

= constant = kn

(A.11)

The constants kl , kg , k and kn can be computed from the inlet boundary conditions. The algebraic expression for upd g is given as, 2 Yl Rg2 (ρ − ρg ) g 9 µl F

AN US

upd g =−

(A.12)

M

Stokes model is applicable for bubbly flows with Reb < 2. The radius of the bubble is chosen such that Stokes model is applicable. In the Stokes regime, friction factor F is set to one. Using the definition of the mixture density from Eq. (A.4), Eq. (A.12) can be rewritten as, upd = kd Rg2 Yl αl (ρl − ρg ) g

(A.13)

PT

ED

where kd = −2g/9µl , a constant which can be calculated. Using Eq. (A.5) in Eq. (A.8), the mixture velocity for section 2 can be written as, Yg2 upd kl g2 u2 = + αl2 Yl2

(A.14)

AC

CE

Using Eq. (A.14) in Eq. (A.9) for section 2 gives, 



upd kl g2  ρg2 αg2  + = kg αl2 Yl2

(A.15)

Using Eq. (A.13) for section 2 in Eq.(A.15) gives, 3 2 kd Rg2 ρg2 (ρl − ρg2 ) αl2 − kd Rg2 ρg2 (ρl − ρg2 ) αl2

+ (ρg2 kl + kg ) αl2 − ρg2 kl = 0

32

(A.16)

ACCEPTED MANUSCRIPT

3 , where k = Using the definition for αg2 = kv ng2 Rg2 v equation can be written as,

4π/3,

3 3 ρg2 kv kl ng2 Rg2 + kg ng2 kv Rg2 − kg





3 6 1 − 2ng2 kv Rg2 + n2g2 kv2 Rg2 = 0(A.17)

CR IP T

5 + ρg2 kv kd (ρl − ρg2 ) Rg2 ng2

the above

Similarly, the gas number density equation can be written as, 2 kl ng2 + kd (ρl − ρg2 ) Rg2 ng2

3 + kn ng2 kv Rg2 − kn = 0



3 6 1 − 2ng2 kv Rg2 + n2g2 kv2 Rg2



(A.18)

M

AN US

In equations (A.16) - (A.18), ρg2 is unknown. ρg2 can be calculated using the ideal gas equation if pressure at section 2 (p2 ) is known. p2 can be computed from the mixture momentum equation. For fixed radius assumption, αg2 can be obtained from Eq. (A.16). However, for variable radius, equations (A.17) and (A.18) have to be solved to obtain Rg2 and ng2 and then calculate αg2 . The mixture momentum equation by neglecting the diffusion in x direction and the contribution of the drift velocity can be written as, ∂ρu ∂x

= −

∂p + ρg ∂x

(A.19)

ED

If the distance between section 1 and 2 is small (∆x ), ρ can be approximated by an average between the two sections. Integration of the above equation and using Eq. (A.10) gives, !

PT

k2 ρ

2

−

k2 ρ

!

1

= p1 − p2 + ρavg g∆x

(A.20)

AC

CE

In the above equation, if section 1 coincides with the inlet boundary, then the inlet pressure, p1 is known and ρ1 can be computed from the inlet boundary conditions using Eq. (A.4). Hence, the above equation can be solved for p2 if ρ2 is known. However, ρ2 is dependent on αg2 by Eq. (A.4). Thus, Eq. (A.16)/(equations (A.17) and (A.18)) and Eq. (A.20) have to be solved iteratively to obtain all the variables at section 2. These variables in turn can be used to compute the values at the next section and hence, marching in this way, solution in the complete pipe can be obtained. The iterative solution procedure is outlined as follows: 1. Guess an initial value for ρ2 . 33

ACCEPTED MANUSCRIPT

CR IP T

Compute p2 from Eq. (A.20). Compute ρg2 using p2 in the ideal gas equation. Solve Eq. (A.16)/(equations (A.17) and (A.18)) to obtain αg2 . Compute the new value of ρ2 from Eq. (A.4). Repeat steps 2-5 until convergence.

AC

CE

PT

ED

M

AN US

2. 3. 4. 5. 6.

34