The phase distribution of upward co-current bubbly flows in a vertical square channel

Abstract

In this work one shows experimental data and numerical results of the void fraction distribution in vertical upward air-water bubbly flows in a square cross-section channel. To measure the void fraction distribution one used a single wire conductive probe. The averaged void fraction ranged from 3.3% to 15%; the liquid and the gas superficial velocities varied from 0.9 m/s to 3.0 m/s and 0.04 m/s to 0.5 m/s, respectively. The experimental results for the void fraction distribution were compared with numerical calculation performed by an Eulerian-Eulerian implementation of the Two-Fluid Model. In this work one performs the turbulence modeling with three approaches: using an algebraic model, the k-epsilon two-phase model and the k-epsilon two-phase two-layer model. Comparisons between the experimental and numerical data revealed, in general, good agreement.

Vertical upward bubbly flow; phase distribution in bubbly flows; bubbly flow in square channel

TECHNICAL PAPERS

The phase distribution of upward co-current bubbly flows in a vertical square channel

A. de Matos; E. S. Rosa; F. A. Franca

Multlab -Dept. of Energy, FEM - Unicamp, 13081-220 Campinas, SP. Brazil, [e-mail]

ABSTRACT

In this work one shows experimental data and numerical results of the void fraction distribution in vertical upward air-water bubbly flows in a square cross-section channel. To measure the void fraction distribution one used a single wire conductive probe. The averaged void fraction ranged from 3.3% to 15%; the liquid and the gas superficial velocities varied from 0.9 m/s to 3.0 m/s and 0.04 m/s to 0.5 m/s, respectively. The experimental results for the void fraction distribution were compared with numerical calculation performed by an Eulerian-Eulerian implementation of the Two-Fluid Model. In this work one performs the turbulence modeling with three approaches: using an algebraic model, the k-e two-phase model and the k-e two-phase two-layer model. Comparisons between the experimental and numerical data revealed, in general, good agreement.

Keywords: Vertical upward bubbly flow, phase distribution in bubbly flows, bubbly flow in square channel

Introduction

It is usual that two-phase bubbly flows occur in channels of non-circular cross section. Equipment having these channels are evaporators used in chemical industries, nuclear light-water reactors, gas-liquid separators used in petroleum plants, to mention a few. The phase distribution on the cross-section of channels influences the transport processes. This explains the vast amount of experimental data, modelling and analysis effort focusing the phase distribution in vertical bubbly flows in round pipes. Contrasting, there are not so many reports presenting measurements and numerical calculation for bubbly flows in non-circular channels. One of the reasons is, probably, due to the experimental difficulties in measuring the void distribution in sections that do not have axial symmetry. Other is the fact that bubbly flows in oblong pipes takes the modelling an its numerical implementation to the limit. The constitutive equations expressing interfacial transfer are not fully established and there is not enough experimental data to support turbulent models that are physically robust. Moreover, the implementation of a multidimensional Two-Fluid Model requires computational effort, power and time, to be solved with an adequate level of convergence and resolution. Hence, there are just a few papers presenting numerical results of bubbly flows in non-circular channels and, moreover, comparing them with experimental data.

Most of the work done on gas-liquid flows in non-circular channels focused on the mechanistic modelling of pattern transitions or approach them with one-dimensional models to calculate pressure drop or mean void fraction. These were the subjects of Mishima et al. (1994), Troniewsky and Ulbrich (1984), Hibiki et al. (1994), Mishima et al. (1993) and Keska and Fernando (1994).

Among the papers presenting measurements of the void fraction distribution in channels having corners one may mention Sadatomi et al. (1982), Shiralkar and Lahey (1972) e Moujaes and Dougall (1987). The former generated a vast data bank on vertical air-water flows in rectangular, isosceles-triangular and annular channels, directed towards the calculation of the pressure gradient. The authors mentioned that the void distribution is one of the most important aspects of two-phase flows in non-circular channels. Their measurements included the cross-sectional distribution of void fraction in bubbly flows, which depicted wall and corner peaking. The later depicted data on void fraction distribution, gas and liquid velocity distribution and bubble size in narrow rectangular channels using both an optical fiber probe and a hot film anemometer. Some of their results for the void fraction, however, were presented in terms of the averaged values over the short and long dimensions of the channels, respectively. This procedure filtered the local distribution data and revealed only a limited aspect of the void distribution across the channel section. Non-filtered void distribution data was depicted only in the longer section center-plane. Similarly to the work of Shiralkar and Lahey (1972), who focused on bubbly flows in eccentric annulus, various authors dealt with more complicated flow geometry, like those taking place in nuclear reactor fuel rod bundles. For the sake of shortness, they will not be mentioned herein.

Lahey et al. (1983) and Lopes de Bertodano et al. [1994a, 1994b) measured and modelled air-water vertical bubbly flows in triangular channels. Using hot film anemometers, single and X-probes, they measured the void fraction and the instantaneous velocity components all over the channel cross-section, and took fine grid measurements on the corners. From these they were able to calculate the mean velocity and the Reynolds stress components, revealing the turbulent structure of the flow. Again, the void fraction distribution depicted wall and corner peaking, with values somewhat higher in the corners.

Hoping to improve the understanding of the phase distribution in non-circular square channels, this work presents, discusses and compares experimental and numerical data on vertical upward bubbly flows taking place in a square cross-section channel. The side of the channel was 0.034 m. The flow was at near atmospheric pressure (0.3 bar to 0.55 bar) and ambient temperature (23ºC). To measure the void fraction distribution one used a single wire conductive probe. The averaged void fraction ranged from 3.3% to 15%, for liquid and gas superficial velocities varying from 0.9 m/s to 3.0 m/s and 0.04 m/s to 0.5 m/s, respectively. The measurements depicted the bubbles concentrating close to the walls and a strong void concentration in the channel corners.

To model the flow one used the Two-Fluid Model. A steady-state isothermal flow with no phase change has been considered. The phase distribution over the channel section came out from a balance between the radial pressure force, the lateral lift, the wall force and the turbulence-induced lateral pressure field. One expects the turbulence playing an important role on the phase distribution over the channel cross section. For this reason, three different turbulence models have been tested: an algebraic model, the standard k-e added with a bubble induced-turbulence term and the k-e two-layer model. The numerical results and its comparison with the experimental data are presented and commented. For convenience, Phoenics, a CFD software that utilizes the finite volume approach, solved the system of equations.

Nomenclatures

Latin

C_pi, C_D, C_vm, C_{li =} pressure, drag, virtual mass and lateral lift coefficients, dimensionless

C_e1, C_e2 = constants in turbulent model

C₁, C_mB, C_m = constants in turbulent model

C_pipa, C_la = constants in lateral lift equation

D = bubble diameter, (m)

E = constant

g = gravitational field (m/s²)

k = Von Karman's constant (dimensionless)

l_m, l_e = scale length, (m)

M_k^D = generalized interfacial force, (N/m³)

M_k^s = generalized interfacial force, (N/m³)

N_g = phase density function, dimensionless

p = pressure, (N/m²)

P = kinetic energy production, (J)

r = radial displacement

Re = Reynolds number, dimensionless

t = Time, (s)

u = axial velocity, (m/s)

V = velocity field, vectorial quantity, (m/s)

We = Weber number, dimensionless

x, y = distances from wall, (m)

x,y,z, = cartesian coordinates, (m)

Greek

a = local void fraction, dimensionless

d_ij = delta de Kroenecker, dimensionless

e = dissipation rate, (m²/s³)

k = turbulent kinetic energy, (m²/s²)

r = density, (kg/m³)

n, n_T = kinematic viscosity, turbulent two-phase kinematic viscosity, (m²/s)

n_SI, n_BI = shear induced, bubble induced kinematic viscosity, (m²/s)

s = surface tension, (N/m)

s_a, s_k, s_e = constants in turbulence model

x = bubble radius, (m)

t = stress tensor, (N/m²)

t^Re = turbulent stress tensor, (N/m²)

q = azimuthal position

Symbols

x = rotacional operator

Ñ = gradient operator

< > = indicates mean value

Superscripts and Subscripts

k=G , k=L, ki , w = subscripts indicating the phases (gas and liquid)

and the property at the interface and at the wall

+ = superscript indicating a non-dimensional quantity

D = indicates generalized force

d = indicates drag force

tp = subscript indicating a two-phase property

The Experimental Apparatus and the Measurement Technique

Figure 1 shows a sketch of the experimental set-up. The air-water bubbly flow took place in a vertical aluminum channel of square cross-section, whose side was 34.1 mm. In the mixing chamber the air entered through an inner porous media of cylindrical shape, positioned along the channel axis. This arrangement produced a fairly uniform bubbly flow just upstream the chamber in the full range of air and water flow rates applied. After the mixing chamber the bubbly flow developed along 73 equivalent diameters before reaching the test section. The test section was made of Plexiglas, exactly matching the cross section of the aluminum channel, allowing flow visualization.

Ordinary tap water and air were the working fluids, flowing at nearly atmospheric pressure, from 0,3 barg to 0,55 barg, and ambient temperature of 23ºC. A centrifugal pump circulated the water. The air came from a compressed air tank maintained at constant pressure. To measure the air and the water flow rates within 2% of overall accuracy, one used a laminar element and a vortex meter, respectively. The range of water and air superficial velocities - referred to the channel hydraulic diameter at the test section conditions - applied in the experiments were j_L = 0.9 m/s to 3.0 m/s and j_G = 0.04 m/s to 0.5 m/s, respectively. The mean void fraction in the test section varied from <a> = 3.3% to 15%. The bubble size was in the range of 2.5 mm to 3.5 mm, measured from digital images generated by a fast (1000 frames/sec) movie camera. As the channel had square cross-section, the mean bubble diameter was taken, sampling a number of bubbles from various frozen images of the channel lateral view, by comparison with the channel side.

A fully description of the conductive probe used to measure the void fraction, including the electronics and the data acquisition can be found in Dias et al. (2000) and Matos et al. (1999).

An inner insulated copper wire of 120 mm diameter, was the cathode in the electrical circuit. The anode was the probe stem, i.e., the stainless steel needle that sustained the inner wire. The wire was cut in such a way that exposed a cross section area having, roughly, 1,1 x 10^-8 m². When immersed in water, the exposed tip of the wire established electrical contact with the probe stem through the water. When the tip was inside a gas bubble the electrical circuit was open. This is the typical on-off signal generated by a classical conductive probe, from which one calculates the local void fraction. The void fraction uncertainty does not follow straightforward from experiments like this, as non-intrusive direct measurements were not performed for comparisons. However, it can be acessed from global measurements of liquid and gas superficial velocities, see Dias et al. for details.

Figure 2 is a scheme of the assembly used to measure the local void fraction over a curved surface enclosed by the channel walls. To perform the measurements the probe was inserted into the channel through a small hole and was driven by a mechanism that allowed radial, r, and rotational positioning. The probe was then positioned at 105 points (measuring stations) over this surface. These measuring stations were then projected onto a channel cross-sectional plane – the measuring reduced plane – where they could be specified by the coordinates (x, y) at the axial location z. This resulted in a detailed void fraction mapping, required for a proper averaging process in a strongly non-uniform distribution. A micrometer with a digital read-out measured the radial displacement, r, with an accuracy of ± 0,02 mm. The azimuthal position, q, could be set within ± 1º. The fact that the actual measurements were not performed over a cross-sectional plane is of minor significance since the flow is developed and the axial distance between the most apart measuring stations is quite short.

The hardware for the data acquisition consisted of National Instruments^® AT-MIO16, a 12 bits ISA bus board and a 486 PC computer. The hardware sampled, acquired and stored the signals generated by the probe at every measuring station, at a frequency of 50 kHz, during the period Dt of 10 seconds. To calculate the local void fraction from this raw signal, one applied a threshold level and converted it to a Dirac-delta function. This function - in other words, the phase density function N_g(x,y,z,t) - indicates the occurrence of gas (a bubble) at the measuring station on time t. Thus, the local void fraction, a(x,y,z), comes from the time averaging of N_g(x,y,z,t),

where Dt is the sampling period. Details of the data processing technique can be found in Dias et al (2000). It is worthwhile to mention that using the above-described procedure the number of digitally acquired points used to calculate the void fraction distribution over the channel section amounted to 4.2 x 10⁷.

The Mathematical Modelling

One used a rather complete formulation of the Two-Fluid Model, Ishii (1975) and Drew and Lahey (1979), for the mathematical representation of gas-liquid bubbly flows inside a vertical channel of square cross-section. A balance between the radial pressure force, the lateral lift, the wall force and the turbulence-induced lateral pressure field governs the phase distribution on the channel section. Interfacial transport in two-phase flows is embodied in the mass, momentum and the interfacial transfer equations. In a steady-state isothermal flow with no phase change, these equations reduce, after some definitions and operations, to Eqs. (2), (3) and (4):

The interfacial momentum transfer equation links the phases:

The subscript k indicates the phase (k=G or k=L, gas or liquid) and the subscript _ki defines proper values of the variables at the interface, like the velocity and the pressure of phase k at the interface. The local void fraction is a, r is the density, V is the velocity, p is the pressure, g is the gravitational field, t_k^Re is the turbulent stress tensor, M_k^D is the generalized interfacial force and M_k^s stands for the force due to the surface tension. The viscous tensor in both phases was considered to vanish when compared with t_k^Re.

In order to solve this system of equations for gas-liquid bubbly flows, one has to constitute it with proper equations for M_k^D, M_k^s, the phasic pressure difference, (p_ki-p_k), t_k^Re and other variables. Besides some questioning, the first three terms are reasonably established. The mathematical representation of the two-phase Reynolds tensor, t^Re, however, is far from being settled. It plays an important role in shaping the multidimensional structure of bubbly flows since there is a coupling between the two-phase flow turbulence and the phase distribution. In a developed flow, for example, the radial force balance accounts for the radial components of forces embedded in the generalized interfacial force, the radial pressure field and the turbulence-induced lateral pressure field. As mentioned by Lee (1988), if the subtle changes in the void distribution are to be modeled, the non-uniform distribution and the anisotropy of liquid phase turbulence have to be taken into account.

The Reynolds Stress Tensor

As the gas is much less dense than the liquid, as is in bubbly flows herein reported, it is enough to constitute the Reynolds, or turbulent, stress tensor in the liquid phase, t_L^Re. Every known model available to constitute it has limitations and drawbacks. Some are not robust enough to represent the physics of the phenomenon, others bring additional complexities to the system of equations, failing the convergence of the numerical solution. If one adds the fact that some disagreement persists in regard to experimental data provided by different authors, it is clear that the subject of turbulence modelling in two-phase bubbly flows in far from being settled. Being aware of this, in this work one accomplished it with three approaches: using an algebraic model, the k-e two-phase model and the k-e two-phase two-layer model. So, the Reynolds stress tensor is written according the Boussinesq eddy-viscosity approximation:

where k is the turbulent kinectic energy and n_T is the two-phase turbulent kinematic viscosity, respectively.

Following Bertodano et all (1994c), one assumes that the two-phase turbulent viscosity is a superposition of the viscosity due to the Reynolds stress, the shear-induced viscosity, n_SI, plus a bubble induced viscosity, n_BI:

The bubble-induced turbulent viscosity stems from the mixing length concept,

where C_mB is equal to 0.69 and d is the mean bubble diameter.

In order to solve the phasic mass and momentum equations one needs to define closure laws for the shear induced viscosity. Three approaches were used to calculate n_T: an algebraic relation, mentioned herein as the l-vel model, the standard k-e model and the k-e two-layer model. The l-vel model is said not being robust, i.e., it does not represent the physics of the phenomenon. Specifically, the signal change of the Reynolds stress over the pipe diameter, originated by the non-centered peak of the velocity profile that experimental data have shown. It is, however, easy to implement and does not impose complexities, like additional differential equations that could lead to systems that do not converge. The k-e is, nowadays, a turbulence reference model, with a widespread use in single-phase flows. However, the standard k-e model is referred to not giving a good representation of the Reynolds stress tensor close to channel walls, i.e, when the Reynolds number is low, besides requiring a very high numerical resolution due to the steep gradient of dissipation rate of turbulent energy, e, in these regions. In order to improve computational efficiency and convergence rates, as well as to introduce established length-scale distribution in the region near the walls, one alternative is to use the standard k-e model in the fully turbulent region and an algebraic expression for e in the damping layer. This is known as the k-e two-layer model, the third approach one used to calculate n_Ttp.

The algebraic model (l-vel).

The algebraic model is set by the non-dimensional relationship between the distance from the wall, y⁺, and velocity, V⁺. The classical Spalding's (1994) law of the wall is:

where the constants k (von Karman) and E are equal to 0.417 and 8.6, respectively, and the non-dimensional variables are

The derivative of y⁺ in regard to V⁺ is the non-dimensional kinematic viscosity, n_SI⁺, which, by its turn, is . Hence, differentiating Eq. (6) in regard to V⁺ gives rise to the turbulent viscosity:

The k- e model.

Now the shear induced viscosity is set by the k-e model,

where k is the turbulent kinetic energy, e is the dissipation rate and the constant Cm is equal to 0.09. The turbulent kinetic energy and the dissipation rate come from the simultaneous solution of the two-phase version of the the standard k-e model [19]:

In these equations, s_k = 1, s_e = 1.314, s_a = 1, Ce₁ = 1.44, Ce₂ = 1.92. P is the kinetic energy production term:

The single phase "law of the wall" equation bridges the fully turbulent flows with the flow near the wall. The values of k and e near the wall come from:

where t_W and y are the wall shear stress and the normal distance from the wall, respectively.

The k- e two-layer model.

The strategy when using the k-e two-layer model, according to Chen (1995), is to distinguish flow regions. The standard k-e model is applied in the fully turbulent region. In the layer near the wall, an one-equation model, the k-e equation combined with an algebraic expression for e , is used:

In the damping layer near the wall the eddy viscosity is

where

The Interfacial Pressure Term

In the liquid phase, the pressure at the interface, p_Li, is

where p_L is the bulk liquid pressure.

This expression results from the solution of a non-viscous flow around a sphere. In this case the coefficient Cp_i is 0.25. Lance and Bataille (1991) suggested values between 0.5~0.7 for actual bubbly flows. Lahey (1990) used C_pi = 1 in distorted bubbly flows. The value used in the present work was a function of void fraction, as suggested by Spalding (1994),

where C_pipa is 0.5. Following Stuhmiller (1977) one considered that the pressure difference is equal for both phases. Using the Laplace equation the relationship between the gas and the liquid pressure at the interface could be stated:

where s is the surface tension and x is the bubble radius.

The Interfacial Forces

The generalized interfacial force acting on the phase k, M_k^D, is usually divided in the drag force, M_k^d, and the non-drag force, M_k^nd. The non-drag force accounts for the lateral forces and the virtual mass force. The drag force arises due to the relative displacement of the bubble. The drag force on a bubble of diameter D is:

According Lahey (1990), the drag coefficient, C_D, depends on the Reynolds number based on the liquid properties. For Antal et al. (1991), conversely, the Reynolds number must be calculated for the gas-liquid mixture properties. In this work one follows the proposition of Kuo and Wallis (1998). According to these authors, the drag coefficient is

where Re is the Reynolds number and We is the Weber number, both expressed as a function of the bubble diameter and the bubble relative velocity, (V_G-V_L). The expressions for C_D are valid for bubbly flows in an infinitum medium. Marié (1987) showed that next to the pipe wall the drag force alters. In this work, however, one does not account for the wall effect on the drag force.

A lateral force, M^L_L=-M^L_G, arises when the bubble displaces through a liquid with a non-uniform velocity distribution. In an upward bubbly flow it acts to displace the bubble toward the wall:

According to Drew and Lahey (1987), the lift coefficient C_L varies between 0.05 and 0.1 in air-water bubbly flows. More recently their co-authors have suggested greater values: Antal et al. (1991) used values in the range 0.01 ~ 0.5 and Wang et al. (1987) adopted 0.5. The correlation used in this work was due to Spalding (1994)]:

One varied C_La in order to find the most appropriate value fitting the experimental data and found C_La = 0.5. The resulting values for C_L were between 0.22 and 0.48.

In bubbly gas-liquid flows the bubble expansion eventually changes the relative displacement between the gas and the liquid, giving rise to a virtual mass force. The virtual mass force is, Drew and Lahey (1987):

As the virtual mass coefficient, C_vm, relates with displaced volumes, one expected it to depend upon the void fraction, a. However, various authors have used constant values for C_vm. Drew and Lahey (1987) suggested 0.5; Kuo and Wallis (1988) used values within the range 2.0 ~ 3.0, and Lance and Bataille (1991) adopted values ranging from 1.2 to 3.4.

Similarly with has been done with the lift coefficient, the virtual mass coefficient has been adjusted to fit numerical and experimental data. The values that have assured the best agreement, i.e., the "ideal values" of C_L and C_vm, will be subject of prospective comments.

The Numerical Method

The computations were performed in a workstation SUN Ultra10, running the CFD PHOENICS 3.2, the package developed by Spalding (1994). PHOENICS is a transport equation solver with two-fluid model capability. It is based on the SIMPLE algorithm. The IPSA algorithm simulated the interfacial momentum exchange.

A computational grid was set up with 20x20x50 nodes, corresponding to the coordinates (x,y,z), respectively. A marching procedure, exploring the parabolic nature of the flow, was used to solve the conservation equations. The initial condition, at the mixer position, matched up to uniform velocity field and void fraction distribution; the pressure was fixed, indirectly, through the setting of the gas density. The calculation then proceeded along 73 equivalent diameters until reaching the measuring reduced plane.

When working with the k-e turbulence model, careful attention was given when setting the first control volume for numerical integration, in order to follow the constraints imposed by the "law of the wall". With the Reynolds number ranging from 30000 to 40000, the outermost face of the first control volume was located at 2.5 mm from the wall. This makes the position of the first measuring station different from the first numerical volume. Hence, the numerical and the experimental results at the measuring station closest to the wall, in every xz plane, should not, ideally, be compared. When working with the algebraic model and the k-e two-layer model this limitation does not exist: the numerical grid could be refined to rather small distances from the wall.

Regarding the numerical convergence, one adjusted the linear under-relaxation parameter, LINRLX, used in PHOENICS. This parameter varies between 0 and 1 and updates the current variable by a percentage of the previous value. If LINRLX is equal to one there is no under-relaxation; if LINRLX is equal to zero, the variables are not updated between two steps of the integration process. Due to the strong coupling between the phases, the set of equations had to be under-relaxed at different ratios. The pressure used a 0,8 factor, the void fraction was under-relaxed at 0.4; the velocities in the principal and secondary directions were under-relaxed at 0.7 and 0.01, respectively. The other flow variables, such as the turbulent kinematic viscosity and the relative velocity were under relaxed at 0.15. These values might not be the optimum under-relaxation factor, but assured a steady and converging solution. To reach this converged solution, every operational condition, established by the gas and liquid superficial velocities, required 5 hours of computation.

Experimental and Numerical Results

Initial tests were conducted to access the verticality of the channel. The void fraction distribution over the channel cross-section, appearing in the contour map in Fig. 3, showed that reasonable vertical positioning had been achieved. The software Surfer^® produced the contour map.

Even using a less refined measuring grid, these first measurements pointed out to results that occurred in every operational condition: bubbles concentrating near the wall, specially on the channel corners. The void peaking near the wall, typical in vertical upward bubbly flows in round pipes, also occurred in the square cross-section channel.

Definite and more detailed measurements were performed on half of the test section. Exploring the axis-symmetrical nature of the flow, the local void fraction was measured in 105 stations for various operational conditions. The measuring grid was particularly refined in the channel region close to the wall. Corresponding numerical results for the void fraction distribution have been computed. The numerical results refer to: (i) the three turbulent models described in a previous section, and (ii) specific formulations for the lift and virtual mass forces with "ideal values" for C_L and C_VM, as discussed.

Figure 4 is a set contour maps showing the experimental and numerical distribution of the void fraction over the channel cross-section. The maps have been plotted using a commercial package, Surfer^®. Due to the fact that one does not have control over the interpolation procedure embedded in the package, it is not possible to state that they are a true representation of the data. However, they are a good overview of the void fraction distribution over the channel section and permit a quick, besides rough, comparison between the experimental and numerical data.

To the left, the outermost map represents the experimental data. The three subsequent maps to the right stand for the numerical results, each one representing a turbulence model, k-e two-layer, standard k-e and algebraic l-vel. Every sub-set stands for a pair of superficial velocities: the upper, #1, (J_L= 2,12 m/s; J_G = 0,13 m/s); in the middle, #2, (J_L = 1,33 m/s; J_G = 0,14 m/s) and the lowest, #3, (J_L= 0,84 m/s; J_G = 0,03 m/s). The respective mean void fractions are: #1: <a> = 5.8% ; #2: <a> = 9.5% and #3: <a> = 3.4%. The numerical results were obtained, after extensive calculations, for "ideal values" of the lift and the virtual mass coefficients, i.e., the ones that fitted the best the void fraction distribution over the channel wall.

Looking at the contour maps one can state that, in a broad sense, the present implementation of the Two-Fluid Model was able to represent the phase distribution of an ascending vertical bubbly flow in a square cross-section channel. Moreover, the k-e two-layer turbulence model, among those used as constitutive equations, gave the best representation of the phase distribution in the range of the applied superficial velocities. The "ideal value" for the lateral lift coefficient depended on the local void fraction and was in the range C_L = 0.22 ~0.48; for the virtual mass coefficient, was C_VM = 2.

The analysis of the Cartesian plots appearing in Fig. 5, 6 and 7 provided "de-facto" comparisons between experimental and numerical data. Every plot in the figures refers to the void fraction distribution in a transversal y-z plane defined by a x position, each figure relating to a pair of (J_L;J_G). The lowest the x value, the closest to the lateral wall the y-z plane is. The highest x value, x= 17.35 mm, defines the mid-transversal plane.

When one looks at the void distribution on the two planes close to the wall, x= 0.8 mm and x = 1.8 mm, it is clear that the numerical distribution, in general, did not have a good agreement with the experimental data. In these planes, ning from corner to corner, the experimental void fractions were, in general, greater than the numerical results. The difference amplified as the mixture velocity decreased, from Fig. 5 to Fig. 7. Also, the experimental distribution showed subtle changes the numerical calculations were not able to pair. The statement holds true no matter the turbulence relation embedded in the Two-Fluid Model and the magnitude of the gas and liquid superficial velocities. There are several factors, different in nature, which could explain the discrepancies between experimental and numerical results.

The experimental distribution lacked symmetry in some cases, pointing out to experimental difficulties in measuring the void fraction distribution at locations very close to the lateral wall or in the channel corners. The probe positioning could be argued as one of these difficulties. Furthermore, at these locations, where the velocities are lower and secondary flows exist, the probe might interfere with the flow field. If one looks at the similitude between the experimentation and the mathematical representation of the flow, new problems arise. The actual distribution of the void fraction at locations close to the wall depends on a "physical" characteristic of the flow, the bubble size distribution function. The Two-Fluid Model does not have the bubble size as a constraint when calculating the void fraction distribution. In other words, in real world the gas volume has a finite dimension – the probe only "sees" bubbles greater than a certain size, for example, while the Two-Fluid Model deals with the mixture as a continuum. Hence, it might not make sense to compare experimental and numerical results in locations close to the wall if the calculation grid was smaller than the bubble size, for example. This was the case of the measuring stations on the first transversal planes, which had x = 0.8 mm and x = 1.8 mm.

When one analyses the experimental and numerical void distribution on inner planes, from x= 3.5 mm to x = 17.35 mm, a different figure appeared. Good agreement was achieved, no matter the magnitude of the mixture velocity or the mean void fraction. Among the turbulence models used, the k-e two-layer gave, consistently, the best results. On inner planes the void fraction profile was quite flat in the channel center region and the numerical results delivered by the Two-Fluid Model with the k-e two-layer model embedded matched up to it. Getting closer to the walls the experimental and numerical profiles were somewhat detached. However, the numerical results delivered by the Two-Fluid Model with the k-e two-layer model embedded had the highest gradient. Thus, if not agreed with the experimental data in shorter distances from the wall, at least had a similar trend.

The numerical results of the Two-Fluid Model with the standard k-e and the algebraic l-vel models embedded did not pair the experimental profiles as close as the solution provided by the k-e two-layer model. The numerical results were higher than the experimental ones in the channel center region, and were lower close to the wall.

Conclusion

A vast number of measurements of the void fraction distribution in an air-water vertical bubbly flow in a square cross section channel has been made and disclosed. To perform the measurements a single wire conductive probe has been used. The experimental data revealed that the void fraction profiles presents the wall peaking that is typical in upward bubbly flows in round pipes. Moreover, the void concentration in the pipe corners was particularly high. The void fraction profile has been calculated using a rather complete implementation of the Two-Fluid Model, constituted with three turbulence models: an algebraic l-vel model, the standard k-e and the k-e two-layer models.

The Two-Fluid Model with the k-e two-layer model embedded gave the best representation of the void fraction distribution. Close to the walls and on the channel corners the experimental and numerical profiles were somewhat detached, in some cases. However, the solution provided by the k-e two-layer model, having the highest void gradient next to the wall, showed the same trend as the experimental data.

Acknowledgements

The authors wish to acknowledge the finnancial support provided by FAPESP, which granted Dr. de Matos a scholarship, and FINEP and CNPq, for the financial aid for instrumenting the multiphase flow loop. They all are sponsoring agencies in Brazil.

References

Antal S.P., Lahey Jr, R.T. e Flaherty, J.E., 1991, Analysis of phase distribution in fully developed laminar bubbly two-phase flow, Int. J. Multiphase Flow, vol. 17, n. 5, pp. 635-652.

Chen, K., 1995, Comparison of different k-e models for indoor air flow computations, Num. Heat Transfer part B, vol. 28, pp. 353-369

Dias S.G., França F.A. e Rosa E.S., 2000, Statistical method to calculate local interfacial variables in two-phase bubbly flows using intrusive crossing probes, Int. Journal of Multiphase Flow, 26, pp. 1797-1830.

Drew, D. and Lahey Jr, R.T., 1979, Application of general constitutive principles to the derivation of multidimensional two-phase flow equations, Int. Journal Multiphase Flow, vol. 5, pp. 243-264.

Drew, D.A. and Lahey Jr, R.T., 1987, The virtual mass and lift force on a sphere in rotating and straining inviscid flow, Int. Journal Multiphase Flow, v.13, n. 1, pp. 113-121.

Hibiqui T., Mishima K., Yoneda K. Fujine S. Tsuruno A. Matsubayashi M., 1994, Visualization of fluid phenomena using a high frame-rate neutron radiography with a steady thermal neutron beam, Nuclear Instruments and Methods in Physics Research a vol. 351 pp. 423-436

Ishii, M., 1975, Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris - France

Keska, J.K. e Fernando, R.D., 1994, Average physical parameters in an air-water two-phase flow in a small, square-sectioned channel, J. Fluids Engineering, Vol. 116, pp. 247.

Kuo, J.T. and Wallis, G.B., 1988, Flow of bubbles through nozzles, Int. J. Multiphase Flow, Vol. 14, n. 5, pp. 547.

Lahey R.T., 1990, The analysis of phase separation and phase distribution phenomena using two-fluid models, Nuclear Engng. & Design, Vol. 122, pp. 17.

Lahey, Jr., R.T., Lopez de Bertodano, M. and Jones, O.C., 1993, Phase distribution in complex geometry conduits, Nuclear Engineering and Design, vol. 141, pp. 177-201.

Lance, M. and Bataille, J., 1991, Turbulence in the liquid phase of a uniform bubbly air-water flow, J. Fluid Mechanics, vol. 222, pp. 95-118.

Lee, S. J., Turbulence modelling bubbly two-phase flows, PhD Thesis, Rensselaer Polytechnic Institute, USA, 1988.

Lopez de Bertodano, M.L., Lahey, R.T. and Jones, O.C., 1994, Phase distribuiton in bubbly two-phase flow in vertical ducts, Int. J. Multiphase Flow, vol. 20, no. 5, pp. 805-818.

Lopez de Bertodano, M.L., Lahey, R.T. and Jones, O.C., 1994, Turbulent bubbly two-phase flow data in a triangular duct, Nuclear Engineering and Design, vol. 146, pp. 43-52.

Lopez de Bertodano, M.L., Lahey, R.T. and Jones, O.C., 1994c, Development of a k-e model for bubbly two-phase flow, Trans. ASME, vol. 116, march, pp. 128-134.

Marié, J.L., 1987, Modelling of the skin friction and heat transfer in turbulent two-component bubbly flows in pipes, Int. J. Multiphase Flow vol. 113, nº 3, pp 309-325

Matos, A., Rosa, E. S., Franca, F. A., Morandim, M. and Morales, R., 1999, Vertical upward bubbly air-water flow in a vertical square section channel (in portuguese), Proceedings of the XV COBEM (CD-ROM), Águas de Lindóia, SP, Brazil.

Mishima, K. Hibiki, T. e Nishihara, H., 1993, Some characteristics of gas-liquid flow in narrow rectangular ducts, Int. J. Multiphase Flow, Vol. 19, # 1, pp. 115-124. Mishima, K., Hibiki, T., Yoneda, K., Fujine, S., Kanda, K., Nishihara, H. e Chung, N., 1994, Annu. Rep. Res. Reactor Inst. Kyoto University, Vol. 27, pp. 245-253.

Moujaes, S e Dougal R. S., 1987, Experimental investigation of cocurrent two-phase flow in a vertical rectangular channel, Canadian J. Chem. Eng. Vol 65, October.

Sadatomi M., Sato Y., e Saruwatari S., 1982, Two-phase flow in vertical noncircular channels, Int. J. Multiphase Flow, vol. 8, n. 6, pp. 641-655.

Shiralkar, B. and Lahey, Jr., R. T.1972, Diabatic local void fraction measurements in Freon 114 with hot-wire anemometer, ANS Transactions, 15(2).

Spalding, B.D., "The Phoenics Encyclopaedia", 1994.

Stuhmiller, J.H., 1977, The influence of interfacial pressure forces on the character of two-phase flow model equations, Int. J. Multiphase Flow, Vol.3, pp 551.

Troniewski, L. e Ulbrich, R., 1984, Two-phase gas-liquid flow in rectangular channels, Chem. Engineering Science, Vol. 39, # 4, pp. 751-765.

Wang S.K., Lee S.J., Jones Jr O. C., e Lahey Jr R. T., 1987, 3-D turbulence structure and phase distribution measurements in bubbly two-phase flow, Int. J. Multiphase Flow, vol. 13, n 3, pp. 327-343.

Paper accepted September, 2004. Technical Editor: Aristeu Silveira Neto.

Publication Dates

• Publication in this collection
  07 Oct 2004
• Date of issue
  Sept 2004

History

• Accepted
  Sept 2004
• Received
  Sept 2004