Modeling of Turbulent Bubbly Pipe Flow with the CFX4.2 Computer Program

 

ANDREY A. TROSHKO and YASSIN A. HASSAN

 

Nuclear Engineering Department

Texas A&M University

College Station, TX, 77832-3133, USA

Phone: 1-(409) 845-7090

Fax: 1-(409) 845-6443

e-mail: y-hassan@tamu.edu

 

 

ABSTRACT

The purpose of this paper is to investigate the ability of the CFX4.2 computer program to model adiabatic, turbulent, upward bubbly flow in a pipe. An effort was performed to estimate the numerical uncertainty of the calculation. Only then would a comparison of the calculations with the experimental results clarify the shortcomings of the physical models. The comparison of the predictions with the experimental data indicates that the code was able to predict effects of liquid velocity flattening and turbulence intensity reduction caused by the presence of the bubbles. It was found that application of Sato's model of bubble induced turbulence led to a better agreement with the experiment. However, the calculation of the wall void fraction peak value was underpredicted.

 

Keywords: bubbly flow, turbulence, two-fluid model, numerical uncertainty, lift force

 

INTRODUCTION

Turbulent bubbly flows are encountered in a wide variety of industrial processes. Their accurate modeling can significantly reduce design costs and improve performance. Methods of computational fluid dynamics (CFD) are one of the tools to achieve such predictions. One of the computational modeling techniques for multiphase flows is the two-fluid model. This model treats phases as interpenetrating media. Distribution of each phase is described by an additional degree of freedom, namely, void fraction. Separate conservation equations are written for this variable in addition to conservation equations for momentum and energy of each phase. The resulting system of interconnected, non-linear equations must be solved. Since most of the multiphase flows are turbulent, an appropriate turbulence model has to be implemented. The development of multiphase turbulence models is currently in its rudimentary state. Recently, several studies have been performed to study turbulence in multiphase flows. Lopez de Bertodano (1992) and Morel and Bestion (1997) reported a successful application of two-phase turbulent models to model bubbly and slug flows in pipes and triangular conduits. Their models are essentially two-phase extensions of the well known and Algebraic Stress models of single-phase turbulence. The CFX4.2 computer code is a commercially available program with extensive capabilities to predict multiphase flows including turbulent flows. The purpose of this study is to verify the code's current mathematical model of turbulent multiphase flow. This is a part of authors' effort to create an advanced model of bubbly turbulent flow. A comparison of the code's predictions against available experimental data (Wang et al. 1987) on upward, adiabatic, turbulent bubbly flow in the pipe was performed. A numerical accuracy study was performed to estimate the numerical uncertainty of the calculations.

 

BRIEF DESCRIPTION OF CFX4.2 CODE

The code predictions are based on the numerical solution of Navier-Stokes equations with various extensions. With respect to multiphase flows, CFX4.2 is able to handle multifluid model when different phases can be incompressible / compressible, turbulent / laminar. It utilizes a non-staggered grid for all variables and a conservative control volume approach to discretize differential equations. Pressure-velocity coupling is handled through SIMPLEC or PISO algorithms. The code is able to calculate steady state and transient flow conditions.

 

DESCRIPTION OF THE EXPERIMENT

Wang et al. (1987) conducted experiments on adiabatic air/water up and down bubbly flows in a pipe (57.15 mm i.d.). Local void fraction, liquid velocity and Reynolds stress components were measured.

 

MATHEMATICAL MODEL

Numerical solution of the following mathematical model was obtained to compare with experiment. The model solves for the following conservation equations:

continuity equation

(1).

momentum equation

(2).

In Equations (1) and (2): is the local void fraction of phase , is the velocity, is the density, is the effective velocity, is the gravity force, is the pressure, is the interfacial force which included drag, lift, turbulent dispersion, wall lubrication and virtual mass forces. Both phases were assumed to be turbulent. The effect of turbulence was described by the model of turbulence with addition of the bubble induced turbulent viscosity model of Sato et al. (1981):

 

(3)

where subscripts refer to liquid and gas phase respectively, is the bubble diameter, are empirical constants. The conservation equations for turbulent kinetic energy and dissipation rate are obtained by summing up respective equations for each phase and coincide with standard single phase form:

 

(4)

where are standard single phase net sources containing production and dissipation terms. In Equation (4) are empirical Prandtl numbers, and

 

COMPUTATIONAL CONDITIONS AND ASSUMPTIONS

Numerical solution of the mathematical model described above was obtained on a 2D computational domain. Due to the axial symmetry only half of the pipe diameter was simulated. Inlet profiles for phase velocities and void fraction were not provided in the experimental reference (Wang et al. 1987). Thus, constant inlet values were assumed. In conservation equations, the order of approximation is determined mainly by the advective term approximation order. In all computations, all variables were approximated by the TVD second order scheme with the Van Leer (1977) flux limiter to ensure boundedness of solution. Linear interpolation was used in the pressure correction equation. Outlet boundary condition was located 35 diameters downstream and imposed zero axial gradients on all variables except pressure. This corresponds to the fully developed flow in accordance with experimental data. Solid wall boundary conditions for velocities and turbulence scales were implemented through standard wall function approach. Zero flux through the solid wall was prescribed for void fraction.

 

NUMERICAL ACCURACY STUDY

In order to establish a meaningful comparison between CFD prediction and experiment, numerical uncertainty should be first estimated. This uncertainty primarily consists of discretization, convergence and outlet boundary condition errors. Round-off error is usually negligible. In all calculations, the convergence criterion was used when total mass error dropped below of total inlet mass flow rate. It was noticed that this corresponded to a mass error for each phase below for the inlet mass flow rate. The value of each variable at the monitoring point (middle of domain) also seized to change. Thus, error associated with convergence does not exceed 1% for each variable for all calculations. To estimate discretization error several calculations with different refinements were performed. Number of control volumes was 300 in axial directions and 10, 20, 30, 40 and 50 radial nodes. Figures 1 and 2 show results of grid refinement for void fraction and gas axial velocity fully developed profiles. As shown, starting from 40 nodes results essentially do not vary. To check sensitivity to axial refinement, the number of axial nodes was increased to 400 and no significant change in results was found. For second order method accuracy, discretization error associated with grid of size h can be estimated as (Celic and Zhang, 1993) , where is the value of interest, is the refinement constant and h is sufficiently small. Calculations showed that this error do not exceed 3% for all variables in fully developed regime. The effect of the outlet boundary error can be checked by imposing a pressure boundary condition as an outlet and comparing the solution obtained with the zero axial gradient outlet condition. Negligible difference was observed. Thus, total computational uncertainty for a given case can be estimated to be about 5% for all variables in fully developed regime. It is necessary to note that all grid refinement calculations were done without Sato model of bubble induced turbulence.

 

COMPARISON WITH EXPERIMENT

Computational grid of 300*50 control volumes was chosen for comparison with experiment. Figure 3 shows fully developed mean liquid velocity profiles for the single-phase case and two-phase up flows. Single-phase data were obtained through Algebraic Stress Model as implemented in the computer program. It is clear that a good agreement is achieved and introduction of Sato model does not exhibit significant difference with respect to the numerical uncertainty (error bars) for liquid velocity. Figure 4 displays fully developed profile for axial diagonal element of Reynolds stress in liquid. The code predictions demonstrated a significant reduction of turbulent intensity, when the bubbly flow was present. This effect can be explained by the fact that presence of bubbles have the influence to make the liquid velocity profile flatter. This leads to a smaller mean liquid velocity gradient and, thus, reduces a source of turbulence. Calculation without Sato model predicted this effect (Figure 3) and also showed reduction of turbulence. However, the reduction was overestimated. Introduction of Sato model essentially enhances turbulence due to the pseudoturbulence caused by wakes of bubbles. Thus, a calculation with Sato model is a compromise between a turbulence reduction due to velocity profile flattening and its enhancement due to wakes of bubbles. Figure 4 shows that introduction of Sato model exhibits change beyond numerical uncertainty and leads to a better agreement with experiment. Figure 5 shows void fraction fully developed profile. Obviously, calculations with and without Sato model underpredicted void wall peak. Wang et al. (1987) showed that for pipe fully developed flow, void fraction distribution is determined by the magnitude of lift force and lateral diagonal element of the Reynolds stress. Empirical coefficient in the lift force may range from 0.01 (low Re number) to 0.5 (inviscid flow). In all present calculations, a value of 0.1 was chosen as recommended by de Bertodano (1992) for this type of flow. However, it seems that the predictions did not produce the right lateral turbulent intensity. There was no measurement of this value for this experiment.

 

CONCLUSION

Comparison between experimental data and CFX4.2 code prediction was conducted after numerical uncertainty has been estimated. Results show that the code was able to predict observed phenomena such as liquid velocity flattening and turbulence reduction caused by bubbles. The code predicted a good quantitative measure of these effects. It was also found that inclusion of the Sato model of bubble caused pseudoturbulence led to a much better agreement with measured axial turbulent intensity. However, it did not influence the mean liquid velocity profile. The calculations underpredicted freestream void fraction and strongly underpredicted the magnitude of wall void peak. This might be attributed to the calculated lateral turbulent intensity. Experimental data is needed to confirm this assumption.

 

REFERENCES

Celik I. and Zhang W. M., "Application of Richardson extrapolation to some simple turbulent flow calculations", Quantification of Uncertainty in Computational Fluid Dynamics, Proceedings of ASME Fluids Engineering Conference, Washington D.C., FED-Vol. 158, 29-38 (1993).

Lopez de Bertodano M. A., "Turbulent bubbly two-phase flow in a triangular duct", Ph. D. thesis, Renselaer Polytechnic Institute, USA, 1992.

Van Leer B., "Towards the ultimate conservative difference scheme, IV. A new approach to numerical convection", J. Comput. Phys., 23, 276-299 (1977).

Morel C., Bestion D., "Study about turbulence modeling in steam water two-phase flows", Proceedings of ASME Fluids Engineering Division summer Meeting FEDSM'97, Vancouver, Canada, 22-26 June, 1997.

Sato Y., Sadatomi M., Sekoguchi K., "Momentum and heat transfer in two-phase bubble flow, I", Int. J. Multiphase Flow, 7, 167-177 (1981).

Wang S. K., Lee S. J., Jones Jr. O. C., Lahey Jr. R. T., "3-D turbulence structure and phase distribution measurements in bubbly two-phase flows", Int. J. Multiphase Flow, 13, 3, 327-343 (1987).

 

 

Figure 1. Grid refinement for fully developed void fraction profile.

 

 

Figure 2. Grid refinement for fully developed gas velocity profile.

 

 

Figure 3. Comparison with experiment for fully developed mean

liquid velocity profile.

 

 

Figure 4. Comparison with experiment for fully developed axial turbulent intensity profile.

 

 

Figure 5. Comparison with experiment for fully developed void fraction profile.