A correlation for interfacial area concentration in high void fraction flows in large diameter channels

Chemical Engineering Science 131 (2015) 172–186

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

A correlation for interfacial area concentration in high void fraction flows in large diameter channels J.P. Schlegel a,n, T. Hibiki b a b

Department of Mining and Nuclear Engineering, Missouri University of Science and Technology, 222 Fulton Hall, 301W 14th St, Rolla, MO 65409, USA School of Nuclear Engineering, Purdue University, 400 Central Dr., West Lafayette, IN 47907-2017, USA

H I G H L I G H T S







An algebraic correlation for interfacial area concentration is derived. Correlation is applicable beyond bubbly flow in large diameter pipes. Includes models for Sauter mean diameter for two groups of bubbles. Data of Schlegel et al. (2012, 2014) is predicted with 22.6% error.

art ic l e i nf o

a b s t r a c t

Article history: Received 17 February 2015 Received in revised form 24 March 2015 Accepted 4 April 2015 Available online 13 April 2015

Two phase flows exist as a part of many industrial processes, including chemical processes, nuclear reactor systems, and heat exchangers. In all of these applications the interfacial area concentration is an important parameter for evaluating the interactions between the phases, including drag forces, heat transfer or chemical reaction rates. Many models for interfacial area concentration exist for dispersed bubbly flows; however this type of flow only exists at relatively low void fractions. Very few correlations exist for the prediction of cap-turbulent, slug, or churn-turbulent flows. In this paper a new correlation for predicting the interfacial area concentration beyond bubbly flows in large diameter pipes is derived using a two-bubble-group method (spherical and distorted bubbles as Group-1 bubbles and cap and churn-turbulent bubbles as Group-2 bubbles) and the two-group interfacial area transport equation. The derivation assumes steady state and fully developed flow, and is based on interfacial area transport source and sink terms for large diameter pipes developed by Smith et al., 2012a. Int. J. Heat Fluid Flow 33, 156–167. The resulting equations can be used to predict the void fraction for each group of bubbles and the Sauter mean diameter for each group of bubbles in addition to the total interfacial area concentration. The model is then benchmarked based on the data collected by Schlegel et al., 2012. Exp. Therm. Fluid Sci. 41, 12–22; Schlegel et al., 2014. Int. J. Heat Fluid Flow 47, 42–56. It is found that the correlation predicts the data for Sauter mean diameter of Group 1 bubbles with RMS error of 23.3% and bias of þ 1.83%. For Group 2 bubbles the RMS error is 24.0% and the bias is þ5.35%. This indicates that the correlation somewhat over-predicts the bubble sizes. In spite of this the prediction error remains reasonable compared to the accuracy of previous correlations, and given that the experimental uncertainty can be as high as 15% for some flow conditions. The RMS error and bias in the total interfacial area concentration are 22.6% and  4.29%, respectively. This is consistent with the overprediction of the Sauter mean diameters, but again is reasonable considering the experimental uncertainty and the prediction error of previous correlations. The model is also able to predict the trends found in the experimental data with varied liquid and gas velocities, representing a large improvement over previous modeling efforts. An expanded database of accurate interfacial area concentration measurements at higher pressures would allow further improvement of the model benchmark and expansion of the range of applicability of the model. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Large diameter Interfacial area Void fraction Churn-turbulent Bubble size

n

Corresponding author. Tel.: þ 1 573 341 7703. E-mail address: [email protected] (J.P. Schlegel).

http://dx.doi.org/10.1016/j.ces.2015.04.004 0009-2509/& 2015 Elsevier Ltd. All rights reserved.

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

1. Introduction Two-phase flows exist as part of industrial processes in many fields. In chemical processing or cooling applications, two-phase flows easily provide a large surface area for efficient chemical reactions or heat transfer. The amount of surface area available for heat or mass transfer, or interfacial geometry, is described quantitatively by the interfacial area concentration. This makes the interfacial area concentration one of the key parameters necessary for predicting the performance of such systems. This importance has been known for decades and has been investigated a great deal. For chemical reaction processes and for heat transfer properties, low void fraction bubbly flow systems have been frequently used (Akita and Yoshida, 1974). For chemical reaction processes this region provides the highest interfacial area concentration, and due to the small bubble sizes in this flow regime the gas phase can react more completely before leaving the system. For heat transfer systems this region allows nucleate boiling, which allows for very high heat transfer rates. Thus most of the work for developing interfacial area concentration models has focused on this region (Hibiki and Ishii, 2001a, 2002; Hibiki et al., 2006; Lin and Hibiki, 2014). However high void fraction churn-turbulent and annular flows are also important in many industrial processes. Phase separation, such as separation of gas and oil in the petrochemical industry, typically occurs in the churn-turbulent flow regime and can be strongly affected by the fluid particle size and shape distribution. Churn-turbulent and annular flows are also commonly found in heat transfer applications such as nuclear reactor systems or industrial boilers. These types of flows have not been investigated or modeled as thoroughly (Hazuku et al., 2007). This is largely due to their complex behavior and the associated modeling challenges. In order to evaluate phase separation or the transfer of mass, momentum and energy it is possible to use a two-fluid or driftflux model coupled to a correlation for interfacial area concentration (Ishii and Hibiki, 2010). For these high-void-fraction processes current modeling methods utilize very simple ad-hoc models (Spore et al., 1993; Thermal Hydraulics Group, 1998). Unfortunately these models are not physically realistic and result in very high error when predicting two-phase flow systems. Recent efforts have been made to develop an interfacial area transport equation for use in various flow regimes (Hibiki and Ishii, 2000a; Fu and Ishii, 2003a,2003b; Sun et al., 2003; Smith et al., 2014) or to use various population balance approaches (Lehr and Mewes, 2001; Yao and Morel, 2004; Sari et al., 2009; Huh et al., 2006; Liao et al., 2011). However these transport equations require detailed a priori knowledge of the inlet conditions, either from another simulation or from experimental data, and such information may not be available for some applications. They are also computationally intensive, requiring complex computational systems and the inclusion of additional transport equations in the computational approach. This does not lend the transport equation approach to some engineering design and analysis situations where such a high degree of accuracy and detail is not necessary. Often simplified models can be used to provide more economical solutions or early in the design process to make informed decisions. To address these needs, a correlation should be developed based on physical arguments to capture the behavior of the interfacial area concentration. Some progress has been made in this area by Ozar et al. (2012), however their model is empirical and may not reflect the physical dependencies of two-phase flows completely. There is also the issue of flow channel geometry. The analyses present in the literature have largely focused on interfacial area concentration in round tubes with small diameter (Hibiki and Ishii, 2001a, 2002; Hibiki et al., 2006). Small diameter tubes have

173

a dimensionless diameter that satisfies the condition (Schlegel et al., 2010) Dh Dnh ¼ qffiffiffiffiffiffiffi r 18 σ

ð1Þ

g Δρ

This corresponds to about 5 cm for air–water flows at atmospheric pressure and temperature. These types of flow systems are characterized by the presence of a stable slug flow regime, with regularly-shaped Taylor slug bubbles. Some work has also been performed for annulus flow channels (Ozar et al., 2012), which are often used to simulate nuclear reactor rod bundles. However many industrial systems have dimensionless diameters which are much larger, ranging from 55 to as much as 500. These systems are considered to be large diameter systems (Brooks et al., 2012). Large diameter systems are characterized by the inability to sustain stable Taylor slug bubbles which occupy the entire crosssectional area of the channel. Almost no work has been done in this area due to the complexity of analyzing such flows, despite their importance in many engineering systems. Dimensionless pipe diameters between 18 and 52 are considered a transition region between small diameter and large diameter systems. The pipe diameter can have a significant effect on the interfacial area concentration for void fractions higher than 0.3. In small diameter pipes the slug bubble relative velocity increases as the pipe diameter increases (Ishii, 1977). This reduces the void fraction, thereby resulting in decreased interfacial area concentration. In the transition region between small and large diameter the slug bubbles begin to become unstable, resulting in a higher bubble number density and smaller cap bubbles. This leads to increased interfacial area concentration. In large diameter systems the diameter is expected to have negligible effect on the interfacial area concentration, as very little change is expected in the bubble behavior at various diameters. For dispersed bubbly flow conditions where no cap or slug bubbles are present the pipe diameter is expected to have very limited influence on the bubble diameter. This effect will be mostly through changes in turbulence behavior as described by the Reynolds number. In order to address these challenges research has been undertaken to accomplish several major tasks: (1) Briefly summarize some unique characteristics of large diameter pipe flows; (2) Analyze and evaluate the performance of existing interfacial area correlations; (3) Develop an interfacial area correlation applicable beyond bubbly flow conditions; (4) Benchmark the correlation using data collected in large diameter channels. 2. Two-phase flows in large diameter channels In large diameter channels, stable Taylor cap/slug bubbles are unable to occupy the entire cross-section of the flow channel. When a cap bubble approaches a certain maximum diameter the surface of the bubble is unable to support the weight of the water above it due to Taylor instability. Generally this maximum bubble size is given as a dimensionless diameter of 30 to 52 (Kataoka and Ishii, 1987; Brooks et al., 2012). If the pipe size is smaller than this value, then the bubble occupies the entire cross-section before reaching the stability limit and can form a stable slug bubble. If the pipe is larger than this limit then any bubbles formed which are larger than the stable size limit will be subject to instability and break apart. It should be noted that it is still possible to form bubbles larger than the stable size limit through coalescence of large bubbles. However the resulting bubble will be unstable and will quickly disintegrate rather than form a stable gas slug.

174

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

This difference may seem small; however it has a very large effect on the behavior of the two-phase mixture. The bubble number density for Taylor cap/slug bubbles is higher in large diameter channels. This can cause enhanced bubble-induced turbulence (Serizawa and Kataoka, 1990; Ohnuki and Akimoto, 2001). The enhanced turbulence, in turn, may cause decreased average bubble size and increased interfacial area concentration in larger channels. Flows in large diameter channels are also much more sensitive to inlet conditions than flows in small diameter channels (Hibiki and Ishii, 2001b, 2003). In small diameter channels the relative velocity of the gas phase is strongly dependent on the channel diameter. This is because the channel diameter determines the slug bubble size and therefore the force due to buoyancy. By contrast, in large diameter channels the relative velocity is nearly constant because the cap/slug bubble size is independent of the channel size (Kataoka and Ishii, 1987). There are also significant differences in the interactions between the gas and liquid phase. This can give rise to a negative or zero liquid velocity, or stagnation region, near the wall of the pipe as documented by Ohnuki and Akimoto (2000). All of these characteristics indicate that large diameter channels can be considered a basic form of two-phase flows, as it has negligible restriction due to the system geometry. This makes flows in large diameter pipes an ideal starting place for the development of models, as the interactions between bubbles are largely restricted only by the physics of interfacial interactions. It can then be postulated that models for more complex geometries can be developed by applying the proper assumptions and boundary conditions to models developed for large diameter channels. For this reason, the current study focuses on interfacial area concentration in large diameter pipes.

3. Interfacial area concentration data for model development A number of data sets are available that include interfacial area transport data for conditions beyond bubbly flows, or for flows with void fraction larger than 0.3 (Lin and Hibiki, 2014). This ranges from flows in small diameter pipes (Fu and Ishii, 2003a,2003b) to flows in large pipes (Smith et al., 2012a; Schlegel et al., 2012, 2014). In addition, measurements have been made in rod bundles (Yang et al., 2013) and annuli (Ozar et al., 2008). The effect of the spacer grids present in rod bundles is not well understood and can have significant implications for

characterization of turbulence. Turbulence is a key parameter in evaluating the interfacial area concentration. In addition, the presence of fuel rods can distort the shape of large bubbles. For these reasons, the rod bundle data is not suitable for evaluating an interfacial area concentration correlation at this time. However there are some similarities between rod bundles and large diameter channels, so the data may be used to evaluate a model for applicability in rod bundles. In the current effort the development of an interfacial area concentration correlation is focused on large diameter channels. In light of this, the data of Fu and Ishii (2003a, 2003b), which was collected in small diameter pipes, and the data of Ozar et al. (2008), collected in an annulus, are not suitable for the current comparison. They may be very useful for further model development in additional geometries. Smith et al. (2012a) did perform experiments in pipes with 0.102 m diameter and 0.152 m diameter, however the 0.102 m diameter channel does not exhibit large diameter channel behavior and data in the 0.152 m diameter channel is largely restricted to bubbly flows. Other research including detailed measurements of two phase flows beyond bubbly flow have been performed by Prasser et al. (2007) and Lucas et al. (2010). In these studies the bubble size distributions were reported, and only for a limited number of conditions. Unfortunately, while these studies are valuable for evaluation of CFD calculations, the area-averaged interfacial area concentration cannot be calculated from the published data. Thus this data cannot be used to evaluate correlations for the area-averaged interfacial area concentration. Based on these discussions, the data of Schlegel et al. (2012, 2014) can be used to evaluate interfacial area concentration correlations for large diameter channels. This represents a concern, as data from only a single group of researchers in two studies does not represent a large database for developing a correlation applicable to systems with varying physical properties. Such systems include high pressure or high temperature systems, as well as systems with varying fluid combinations. However these databases were collected very systematically to cover a wide range of void fractions and flow rates including cap bubbly, cap turbulent and churn-turbulent flows. To improve the evaluation of the models, any correlations should be revisited as additional data on interfacial area concentration beyond bubbly flow in large diameter channels becomes available. The experimental facility used to conduct these experiments is shown in Fig. 1. Experiments were performed at pressures of up to

Fig. 1. Schematic of test facility (Schlegel et al., 2014).

10 1

Flow Regime Transitions (Schlegel et al, 2009) DH = 0.152 m

Superficial Liquid Velocity, jf [m/s]

Superficial Liquid Velocity, jf [m/s]

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

DH = 0.203 m

10 0 Bubbly Annular

Churn-Turbulent Cap-Turbulent

10

-1 -2

10

-1

0

1

2

10 10 10 10 Superficial Gas Velocity, jg [m/s]

10

175

1

Cap-Bubbly Flow Injection Flow Regime Transitions D =0.152 m D =0.203 m D =0.304 m

10

0

Bubbly

Annular

Cap-Turbulent Churn-Turbulent

10

-1

10

-2

-1

0

1

10 10 10 10 Superficial Gas Velocity, jg [m/s]

2

Fig. 2. Test conditions (Schlegel et al., 2012, 2014).

300 kPa and for a variety of gas and liquid superficial velocities, as shown in Fig. 2 (Schlegel et al., 2012, 2014). Gas injection is accomplished using an injector unit, also shown in Fig. 1, which consists of seven porous metal spargers (Schlegel et al., 2012). Each sparger is enclosed in a stainless steel annulus which allows rough control of the average inlet bubble size. For very high gas flow rates, these porous metal elements are replaced by tubes with larger, 3/4″ diameter holes for injection as cap-bubbly flow (Schlegel et al., 2014). The facility includes test sections with three different diameters: 0.152 m, 0.203 m and 0.304 m. Experimental measurements were made using electrical conductivity probes. The electrical conductivity probe was first proposed by Neal and Bankoff (1963) for local void fraction measurements. Kataoka et al. (1986) developed a four-sensor design capable of measuring interfacial area concentration for bubbles with arbitrary shape and deforming interface. The four-sensor probe was later miniaturized by Kim et al. (2000), resulting in the probe diagrammed schematically in Fig. 3. This electrical conductivity probe is capable of measuring interfacial area concentration for bubbles as small as 1 mm in diameter. The signal processing scheme uses the position information of the individual sensors, as well as the phase signal from each sensor, to calculate the interfacial velocity of each interface measured by the sensors. This interfacial velocity is then used to calculate the local, time averaged interfacial area concentration (Ishii and Hibiki, 2010). Bubbles which do not make contact with all four sensors are accounted for using a mathematical transformation based on the average bubble properties. The bubble chord length can be used to categorize the counted bubbles according to size. Benchmark experiments performed by Kim et al. (2000) and Manera et al. (2009) compared the measurements of electrical conductivity probes to detailed image processing results and wire-mesh sensor measurements. These studies found the results to be comparable, with differences of 5% to 10% in the measured values. The uncertainty of the electrical conductivity probe depends on the statistical frequency of bubble contacts, the sensor spacing, the measurement frequency, and the average interfacial velocity of the measured bubbles. The possibility of interface distortion due to the intrusive measurement must also be accounted for. Schlegel et al. (2012, 2014) report a measurement frequency of 30 kHz, with the resulting uncertainty in the interfacial area concentration ranging from 5% at very low interfacial velocities to as much as 25% for interfacial velocities of 15 m/s or more. The void fraction uncertainty is reported to be 5%. The probes were positioned within the flow using a computer-controlled positioning system with an accuracy of 0.01 mm.

Fig. 3. Schematic of four-sensor probe measurement (Kim et al., 2000).

4. Existing interfacial area correlations During development of the two-fluid model several semiempirical correlations were developed to estimate the interfacial area concentration. This led to the development of analyses based on two bubble groups, as illustrated in Fig. 4. In these analyses Group 1 is defined as small spherical and slightly distorted bubbles. Distorted bubbles are roughly spherical, but their shape continuously ‘wobbles’ and can become somewhat elliptical or cap-shaped. Group 2 is defined as larger Taylor cap/slug bubbles and churn-turbulent bubbles. The boundary between the groups is determined by the changing drag behavior of bubbles with size, and is defined as the largest possible distorted bubble size (Ishii and Zuber, 1979). The overall modeling strategy for the development of a new correlation is based on the strategy proposed by Ishii and Mishima, 1984. However this equation is modified such that ai ¼ ai1 þ ai2 ¼

6 1α C α  αgs 6α1 6α 2 αgs þ ¼ þ DSm1 1  αgs Dh 1  αgs DSm1 DSm2

ð2Þ

here, and throughout the paper, all quantities are area-averaged. The averaging operators have been dropped for conciseness and clarity. The relations between the total void fraction, group void

176

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

Fig. 4. Illustration of two bubble group analysis for two phase flows (Lin and Hibiki, 2014).

here the void fraction at the bubbly-slug transition and void fraction at the slug-annular transition are given as

αbs ¼

8 > > > > > <





Dh αL ¼ max 0:25  min 1; 22:22



1=2 

g Δρ

σ



; 0:001

G o2000 kg=ðm2 sÞ

2000 αL þ ð0:5 αL ÞG 1000 2000 o G o 3000 kg=ðm2 sÞ > > > > > : 0:5 G 4 3000 kg=ðm2 sÞ

ð6Þ and

0

0 11 !1=2 !1=2 α gDΔρ 3:2α σ g Δρ @ @ αsa ¼ max 0:5; min ; ; 0:9AA ρg jg jg ρ2g ð7Þ

Fig. 5. Illustration of the void fraction in the liquid slug (Ozar et al., 2008).

fractions, and liquid slug void fractions can be summarized as

α ¼ α1 þ α2 ; α1 ¼ αgs

1α ; 1  αgs

α2 ¼

α  αgs α1 ; αgs  1  αgs 1  α2

ð3Þ

The individual group void fractions are correlated using the void fraction in the liquid slug. The liquid-slug void fraction concept is illustrated in Fig. 5. It is essentially a measure of the packing density of small, spherical bubbles in the liquid phase, with the fraction of large cap, slug, or churn-turbulent bubbles neglected. Some of the most well-developed models for two-phase flows have been applied to the prediction of nuclear reactor safety in advanced thermal hydraulic analysis codes. Two of the most current, most commonly used of these codes are RELAP5 and TRAC, reviewed previously by Ozar et al. (2012). Each of these codes utilizes its own method for estimating the interfacial area concentration beyond bubbly flow cases. In the case of RELAP5 (Thermal Hydraulics Group, 1998), the constant C in Eq. (2) is considered to be 4.5. Additionally the first term, which accounts for the effect of small, spherical bubbles, is slightly modified as 3.6αgs/DSm1. The average size of small, spherical bubbles is given by a Weber number criterion as ! DSm1 ¼ 10

σ ρf v2r

ð4Þ

The void fraction in the liquid slug is correlated based on the flow regime transitions and the total void fraction as   α  αbs ð5Þ αgs ¼ αbs exp 8 αsa  αbs

TRAC (Spore et al., 1993) uses a similar model, however the constant C in Eq. (2) is given as 1. Additionally the diameter of spherical bubbles is given as rffiffiffiffiffiffiffiffiffiffi σ DSm1 ¼ 2 ð8Þ g Δρ And the value of the liquid slug void fraction is given by

8 0:3 G o 2000 kg= m2 s > > > <

G  2000 2000 o G o 2700 kg= m2 s αgs ¼ 0:3 þ 0:2 700 > > > 0:5 G 4 2700 kg= m2 s

:

ð9Þ

Ozar et al. (2012) developed another method to calculate interfacial area concentration beyond bubbly flow. Hibiki and Ishii’s (2002) model for the Sauter mean diameter of Group 1 bubbles was used, which is given by    0:335  1=3 4=3   0:239 DSm1 Lo ε Lo ¼ 1:99 : Dh Lo νf

ð10Þ

Then the Group 1 void fraction is given by 8 α α r α1; max > >   <

α1; max  α1;base α1 ¼ α1; max þ α1; max  αcrit α  α1; max α1; max r α r αcrit : > > :α αcrit r α 1;base

ð11Þ The various void fraction limits are given by 8 < 0:235 þ 0:011jfþ jfþ r 6:1 α1; max ¼ : 0:325  0:004jfþ jfþ Z 6:1

ð12Þ

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

177

correlation for the Group 1 interfacial area concentration may also be necessary.

5. New correlation for interfacial area concentration 5.1. Modeling strategy

Fig. 6. Comparison of interfacial area correlations and data in large diameter channels.

αcrit ¼

8 < 0:511 þ0:006jfþ

jfþ r 6:1

: 0:645 0:015jfþ

jfþ Z 6:1

α1;base ¼

8 < 0:099  0:009jfþ

jfþ r6:1

: 0:05

jfþ Z6:1

ð13Þ

ð14Þ

Based on Eq. (2) the individual contributions to the interfacial area concentration can be determined by developing four models: (1) the total void fraction, (2) the Group 1 Sauter mean diameter, (3) the Group 2 Sauter mean diameter, and (4) the void fraction in the liquid slug. The total void fraction is easily determined as part of the solution in advanced thermal-hydraulic codes. For less detailed treatments it can be estimated using a drift-flux correlation (Hibiki and Ishii, 2003). Such correlations are available for many system geometries. One of the most common drift-flux models for small diameter pipes is that of Ishii (1977). Many other exist for various flow regimes and system geometries. Schlegel et al. (2013) recently proposed a set of drift-flux models based on that of Ishii (1977) and Kataoka and Ishii (1987) for large diameter channels. The model proposed by Schlegel et al. (2013) is used to calculate the total void fraction in the current study. The remainder of this section is devoted to the development of models for the remaining three parameters necessary to complete the model. 5.2. Key interfacial area sources and sinks

The Group 2 interfacial area concentration is given using the Ishii and Mishima (1984) model with C ¼1. The predictions given by each of these three models are shown in Fig. 6, which also includes the data collected by Schlegel et al. (2012). The figure clearly shows that the model of Ozar et al. (2012) performs well in bubbly flows, where it is largely reduced to the Hibiki and Ishii (2002) model. The model also performs reasonably for small diameter channels, as evidenced in that paper. However beyond bubbly flows in large diameter channels the model shows similar deficiencies to those noted in the TRAC and RELAP models. This is expected, as large diameter channels are expected to have higher interfacial area concentration than small diameter channels. This is due to the instability of cap and churnturbulent bubbles in large diameter channels. The instability of large Group 2 bubbles has a number of consequences that contribute to the inability of models for small diameter pipes to predict the interfacial area concentration. First is the inability of those models to account for the elevated Group 2 bubble number density and reduced Group 2 bubble size. These factors act to increase the Group 2 interfacial area concentration. Also, the lack of stable slug bubbles tends to result in a much more gradual transition from dispersed bubbly flow to annular flow with a thin liquid film. While the Group 2 interfacial area concentration is small compared to the Group 1 interfacial area concentration, this contribution will not be negligible. The Group 1 interfacial area concentration also may not be well-predicted beyond bubbly flows. These models all use correlations developed for bubbly flows, where the major factors determining bubble size are random collision and turbulent breakup of Group 1 bubbles alone (Hibiki and Ishii, 2000b). Beyond bubbly flow the Group 1 bubble size is largely determined by the size of the bubbles sheared off around the base of cap bubbles and the coalescence rate of the Group 1 bubbles. Because the sheared-off bubbles tend to be much smaller than the average Group 1 bubble size in bubbly flow, this will result in a decrease in the average Group 1 bubble size beyond bubbly flow and an increase in the Group 1 interfacial area concentration. This means that a new

There are a number of key bubble coalescence and breakup mechanisms in two-phase flows. These categories have been used by researchers such as Wu et al. (1998) and Fu and Ishii (2003a, 2003b) and Smith et al. (2012b) to model interfacial area transport. They are illustrated in Fig. 7. First, random collision is defined as the interaction of two bubbles driven by turbulent fluctuations in the liquid phase. The coalescence rate for this mechanism depends strongly on the number concentration of bubbles and the amount of turbulent velocity fluctuation in the flow. When two bubbles do interact, the collision results in drainage of the liquid film between the bubbles. If the film drains sufficiently, then the two bubbles merge into a single larger bubble. Wake entrainment occurs when one bubble is accelerated into the wake of another, larger bubble. The rise velocity is determined by the relative velocity between the phases and the velocity of the liquid phase. Behind a large, rising bubble the liquid phase is accelerated in a wake region that extends a few bubble diameters behind the large bubble. If another bubble is caught in this wake region, the rise velocity of the trailing bubble is greater than the rise velocity of the leading bubble due to this accelerated liquid velocity in the wake region. This difference in rise velocity means that the trailing bubble will catch up to the leading bubble, resulting in bubble collision and coalescence. Turbulent impact represents bubble disintegration due to interactions of bubbles with turbulent eddies in the liquid phase. Large turbulent eddies simply tend to carry bubbles along without major interaction, while very small eddies do not carry much energy and result in very little effect on the bubble surface. Eddies with length scales similar to the bubble size can readily interact with the interface, however. The energy carried by these eddies is absorbed by bubbles as deformations in the interface. This mechanism competes with production of turbulence due to the relative velocity between the bubble and liquid, and can actually result in suppressed liquid-phase turbulence (Kataoka and Serizawa, 1990). If the turbulent eddy has enough energy, the interface is unable to absorb all of it

178

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

Fig. 7. Two-group source and sink term mechanisms by (a) random collision, (b) Wake entrainment, (c) turbulent impact, (d) shearing-off, (e) surface instability (Smith et al., 2012b).

and the extra energy results in the breakup of the bubble into two separate, smaller bubbles. Shearing-off accounts for bubbles which are pulled off of the base of a cap bubble due to the relative velocity between the cap bubble and the surrounding liquid. This phenomenon has been observed experimentally, but is still not well understood. Various mechanisms have been proposed for modeling shearing-off. The simplest models account for a force balance between the shear forces acting on the bubble and the surface tension force of the interface. More complex models have been proposed based on inviscid Kelvin–Helmholtz instability on the side surface of the cap bubble and gas jet breakup dynamics to predict the frequency and size of the sheared-off bubbles. Whichever mechanism is used to develop a model, the result is the production of many small bubbles from a single large cap bubble. Surface instability accounts for the Rayleigh–Taylor instability in the upper surface of the bubble. This instability defines the maximum length scale for which a heavier fluid can rest upon a lighter fluid in a stable fashion. It depends on the radius of curvature of the bubbles, the surface tension between the fluids, and the density difference between the phases. In modeling this effect it is often assumed that any bubble which exceeds the maximum bubble size dictated by Rayleigh–Taylor instability undergoes instantaneous binary breakup. These bubbles can be created from various coalescence mechanisms or through expansion of a compressible gas phase. This assumption is reasonable for relatively low void fraction flows, when coalescences of large bubbles are relatively infrequent. However it should be noted that at high void fractions the time scale for coalescence of large bubbles may be similar in magnitude to the time scale required for instability to initiate and grow on the bubble surface. In this case, the instantaneous breakup assumption may no longer be valid. This may result in an average measured bubble size larger than the Rayleigh– Taylor limit, as the residence time of bubbles larger than this limit will no longer be negligible. This phenomenon results in rapid creation and destruction of large cap bubbles, which is suspected to cause the onset of churn-turbulent flow in large diameter channels. Along with entrainment of the liquid phase, it also contributes to the transition to annular flow in large diameter channels (Schlegel et al., 2009, 2013). In modeling the bubble size for a two-group flow, detailed in the next few sections of the paper, it has generally been assumed

that the inter-group exchanges of void fraction and interfacial area concentration are the dominant factors in determining the bubble size. These mechanisms include: (1) (2) (3) (4)

Random Collision of one Group 1 and one Group 2 bubble. Entrainment of one Group 1 bubble by one Group 2 bubble Turbulent impact breakup of a Group 2 bubble Shearing-off of Group 1 bubbles around the base of a Group 2 bubble

It is also possible for two Group 1 bubbles to coalesce into a Group 2 bubble through either Random Collision or Wake Entrainment. These mechanisms have been generally neglected here. While they are extremely important in predicting the transition from pure bubbly to cap-bubbly flow, once cap-bubbly flow has been established they tend to be much smaller in magnitude than the four phenomena listed above. The resulting balance is discussed in the next two sections.

5.3. Group 1 Sauter mean diameter Hibiki and Ishii (2002) developed a model for the Sauter mean diameter of bubbles in dispersed bubbly flow, for void fractions less than 0.3, by assuming steady-state conditions. They then simplified the one-group interfacial area transport equation by assuming that coalescence and breakup occur only by random collision and turbulent impact, solving the simplified equation for the interfacial area concentration. This resulted in a correlation for the interfacial area concentration in bubbly flow, which can be used to provide a correlation for the Sauter mean diameter. Unfortunately, this model does not produce accurate predictions for higher void fraction flows where larger cap bubbles are present. In this case the production of Group 1 bubbles by shearing-off becomes a very significant contributor to reducing the bubble size. Their approach can be used however. If a two-group interfacial area transport equation is considered (Smith et al., 2012a), then for Group 1 the steady-state interfacial area concentration balance (neglecting

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

results in:

expansion due to pressure changes) becomes

ϕ

1 RC þ

ϕ

1;2 RC þ

ϕ

1 WE þ

ϕ

1;2 WE þ

ϕ

1 TI þ

ϕ

1;2 TI þ

ϕ

2;1 SO

179

¼0

ð15Þ

N DSm1 ¼

h  1=3 i pffiffi Þ7=4 8 α1 0:151C 1 α2=3 ð1  α2 Þ2=3 þ 2ð1  α11=2 1 1 α2 N vr þ 0:225α1 3 DSm;1 h  i3ð1  α2 Þ ¼ C2  2=5 ρf DSm;2 0:0447C ð1  αÞN vr þ 0:0258N 1

The form of each of these seven terms varies widely, so that this equation cannot be solved explicitly for the interfacial area concentration. It can however be simplified a great deal by linearizing the equation, eliminating higher-order terms. This results in 1;2 2;1 2;1 ϕ1;2 RC þ ϕWE þ ϕTI þ ϕSO ¼ 0

ð16Þ

Solving this equation is a much simpler exercise. Based on the Interfacial Area Transport Equation (IATE) for large diameter channels developed by Smith et al., (2012a) these terms can be expressed as 22:6C RC ε

1=3

α

5=3 1

α

2 4=3 2 DSm;2 DSm;1 þ 11:88C WE

  1=3 vr2 C D2 þ vr1  vr2 α1 α2 DSm;2 DSm;1

¼ 122C TI ε1=3 ð1  αÞDSm;2 D2Sm;1 þ 288C SO 4=3

2=5 ρf3=5 v1=5 r2 σ 3=5 ρg D2=5 N Wec h

α2 D2Sm;1

ð17Þ

This equation can be further simplified. If one defines a Weber number such that

ρf v2r2 Dh N We ¼ σ

ð18Þ

N vr ¼

ε1=3 Lo1=3 vr2

N DSm2 ¼

;

ð19Þ

ρg

We

ð22Þ The constants which have been added to Eq. (22) in this case absorb several of the assumptions made in deriving the equation. 5.4. Group 2 Sauter mean diameter Using the same approach for Group 2 bubbles as described above, the interfacial area concentration balance begins as 1;2 2 11;2 1;2 2 2 2;1 ϕ11;2 RC þ ϕRC þ ϕRC þ ϕWE þ ϕWE þ ϕWE þ ϕTI þ ϕSO ¼ 0

ð23Þ

Substituting the bubble coalescence and breakup kernels developed by Smith et al. (2012a), multiplying by D3Sm2 to eliminate negative exponents on the Group 2 Sauter mean diameter, then linearizing the equation by removing all terms with order greater than 1 results in 1;2 2 2;1 ϕ1;2 RC þ ϕWE þ ϕTI þ ϕSO ¼ 0

or 35:7C RC ε1=3 α1 α22 DSm2 DSm1 þ 33:2C WE ðvr2 C D2 5=3

1=3

1=3

þ vr1  vr2 Þα1 α2 DSm2 þ 7:49C TI ε1=3 ð1  αÞα2 DSm2 DSm1 ¼ 77:8C SO 1=3

and the ratio between the turbulent velocity fluctuation and the relative velocity as

α2

σ

ρf vr2

!

α2

ð24Þ

Then one can define the Weber number as in Eq. (18), the relative velocity number as in Eq. (19), and substitute the various values defined as in Eq. (21). Eq. (24) then becomes

1 DSm2 77:8C SO N We  1=3  o: pffiffi ¼n Þ7=4 Dh ½35:7C RC α1 α2 þ 7:49C TI ð1  αÞN vr2 þ 33:2C WE α1 83 ð1  α2 Þ2=3 þ 2ð1  α11=2 1

ð25Þ

3ð1  α2 Þ

after some algebraic manipulation, Eq. (18) can be rewritten as h  i 5=3 1=3 22:6C RC α1 α22 C 1 Nvr þ 11:8C WE α1 α2 C D2 þ vvr1 1 DSm;1 r2      : ¼ N DSm1 ¼ DSm;2 ρ 122C TI ð1  αÞC 1 Nvr þ 288C SO 3=5α2 2=5 ρf N Wec N We

g

ð20Þ Eq. (20) shows that the dimensionless Group 1 Sauter mean diameter depends on several dimensionless groupings: the ratio between the Group 1 and Group 2 relative velocity, the ratio between the turbulent velocity fluctuation and the relative velocity for Group 2, the Weber number, and the density ratio between the phases. Substituting the values of the constants given by Smith et al. (2012a) and specifying the drag coefficients and relative velocities of the bubbles as (Brooks et al., 2012) C D2 ¼ 83ð1  a2 Þ2

vr1 ¼

!1=4 pffiffiffi σ g Δρ 2 ð1  α1 Þ7=4 : 2

ρf

!1=4

vr2 ¼ 3:0

σ gΔρ ρ2f

ð1  α2 Þ1=2

This equation can then be simplified by substituting the constants given by Smith et al. (2012a) such that ( DSm2 C3 ¼ 0:238N vr2 ½1 þ α1 α2  α NDSm2 ¼ Dh N We

ð21Þ

þ4:95α1

!)  1 pffiffiffi  1=3 8 2ð1  α1 Þ7=4 ð1  α2 Þ2=3 þ  1 : 3 3ð1  α2 Þ1=2 ð26Þ

5.5. The liquid-slug void fraction In large diameter channels, the void fraction in the liquid slug is a misnomer. As stable slug bubbles cannot be maintained, there is no liquid slug between the slug bubbles. However the concentration of Group 1 bubbles in the liquid region, excluding the cap bubbles and churn-turbulent bubbles, is still an important parameter for characterizing the number density of Group 1 bubbles and, therefore, the Group 1 void fraction. While it is a misnomer, the term has been used in the past for small diameter channels and the mathematical definition is the same in large diameter channels as in small diameter channels. For this reason the term continues to be used throughout this paper to provide consistency and clarity, although it should be understood that the phrase is not a physically accurate description.

180

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

When modeling the void fraction in the liquid slug, it is important to first identify the appropriate parameters which most strongly influence the flow. In the case of the void fraction group distribution the most important parameters are the turbulent kinetic energy, the void fraction and the fluid properties such as surface tension and density difference. Generally the turbulent kinetic energy is difficult to model, but the effect of turbulence can be estimated using the turbulent dissipation. Hibiki and Ishii (2001a) correlated the bubbly flow Sauter mean diameter using a non-dimensional number that can be described as a turbulence based Reynolds number,  4=3 1=3  Lo ε N Reε ¼ ð27Þ

νf

It then becomes necessary to predict the turbulent dissipation, which is no small challenge in two-phase flows. Generally in simplified models, dissipation is considered to be equal to production of turbulence. In this case the turbulent dissipation can be divided into shear-induced and bubbleinduced components such that

ε ¼ εsh þ εBI :

ð28Þ

The shear-induced turbulence in two-phase flows can be given by (Hibiki and Ishii (2001a))  j dP εsh ¼ ð29Þ ρm dz f riction The frictional pressure drop is then given by  dP ρ v2 ¼ f tp m m dz f riction 2Dh

ð30Þ

Based on the definitions of the mixture terms, using the Blasius equation for the friction factor, and using the mixture viscosity given by Ishii and Zuber (1979), the shear-induced turbulence can then be given by

εsh ¼

j

ρm

f tp

ρm v2m 2Dh

¼

0:316 v3m 2Dh Re0:25 m

ð31Þ

where N Rem ¼

ρm vm Dh ρm vm Dh 1 ¼ ; μm μf 1  α

ρm ¼ ρg α þ ρf ð1  αÞ;

ð32Þ ð33Þ

and vm ¼ ρ

g jg

þ ρf jf ρm

using this data and a hyperbolic tangent function. The result of the curve fit gives the equation (

αgs ¼ min

α

0:63 tanhð0:00145N Reε αÞ

ð36Þ

Although this model was benchmarked with data from large diameter channels, it may be applicable across a wide range of flow geometries provided that turbulent dissipation is properly modeled. It should also be noted that while the value 0.63 was determined from the least-squares curve fit, it may also represent the volume fraction of tightly packed spheres in a cubic lattice. 5.6. Summary of model equations The equations used to complete the solution set are summarized here. These equations are valid for void fractions higher than 0.3. Based on the drift-flux model proposed by Schlegel et al. (2013), qffiffiffiffiffiffiffiffiffiffiffiffiffi 8 > α r0:51 < 1:2  0:2 ρg =ρf qffiffiffiffiffiffiffiffiffiffiffiffiffi ð37Þ C0 ¼

α > : 1 þ 0:2 ρg =ρf 10:49 α Z0:51 This model uses the drift velocity proposed by Kataoka and Ishii (1987), !1=4 !  0:157 ρg σ gΔρ vgj ¼ 0:030 N μf 0:562 ð38Þ 2

ρf

ρf

which is valid for viscosity numbers smaller than 0.0255 and for void fractions higher than 0.3. For a water–air mixture, the value of the viscosity number is 0.0225. Here the viscosity number is defined as N μf ¼ 

ρf σ

μf

qffiffiffiffiffiffiffi1=2

ð39Þ

σ g Δρ

Then the void fraction is calculated from the drift-flux general expression such that.

α¼

jg   : C 0 jf þ jg þ vgj

ð40Þ

These Eqs. (37)–(40) represent closure relations for calculating the total void fraction, and are not part of the newly developed model. Other closure relations are available; however, Schlegel

ð34Þ

It then remains to determine the bubble-induced component of the dissipation. A simple method for calculating the bubbleinduced turbulence was introduced by Hibiki and Ishii (2001a) such that

εBI  jg g

ð35Þ

This has the advantage of being simple and easy to calculate from the flow condition information. This approach has been used in the current work. For more detailed applications it is possible to use detailed formulations based on the interfacial drag. This turbulence formulation has been used to calculate the dissipation Reynolds number for the data collected by Schlegel et al. (2012, 2014). The dissipation Reynolds number and the total void fraction have been used to correlate the void fraction in the liquid slug, and the results are shown in Fig. 8. Hibiki and Ishii (2001a) also included a diameter-dependent term in their model, however, in large diameter channels the effect of the channel diameter is negligible. A least-squares curve fit was performed

Fig. 8. Comparison of liquid slug void fraction model with experimental data.

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

et al. (2013) evaluated several drift-flux correlations and found that these equations provided the best prediction accuracy. Once the total void fraction is calculated, then the dissipation Reynolds number is calculated and used to determine the void fraction in the liquid slug, (

αgs ¼ min

α

ð360 Þ

0:63 tanhð0:00145NReε αÞ

marked model with the experimental data is shown in the next several figures. For clarity, most of the figures shown include only the data collected by Schlegel et al. (2014) in a 0.152 m diameter tube. Fig. 9 shows the comparison of predicted and measured bubble size for both groups of bubbles. There is a significant amount of random error in the prediction, with calculated RMS errors of 23.3% for Group 1 and 24.0% for Group 2. The calculated biases,

and this is used to calculate the individual group void fractions as

α1 ¼ α

1α gs α  αgs

and

α2 ¼ α1  ααgsgs

0

ð3 Þ

Then the Group 2 Sauter mean diameter can be calculated, ( DSm2 C3 ¼ N DSm2 ¼ 0:238Nvr ½1 þ α1 α2  α Dh N We !)  1 pffiffiffi  1=3 8 2ð1  α1 Þ7=4 2=3 ð1  α 2 Þ þ 1 þ 4:95α1 3 3ð1  α2 Þ1=2 ð260 Þ and the Group 1 Sauter mean diameter,

181

εb ¼

 N  DSm;pred;i DSm;meas;i 1X Ni¼1 DSm;meas;i

ð42Þ

for Group 1 and Group 2 are þ 1.83% and þ 5.35%, respectively. This indicates that the correlations do tend to slightly over-predict the bubble diameters. However the phenomena being described are very complex. Generally errors less than 30% are considered reasonable for interfacial area concentration models. Furthermore the experimental data can be categorized into two groups: experiments where the gas is injected as bubbly flow, and experiments where the bubbles are injected as cap-bubbly flow. These two categories result in slightly different behavior of Group

h  1=3 i pffiffi Þ7=4 α1 0:151C 1 α12=3 α2 N vr þ 0:225α1 83 ð1  α2 Þ2=3 þ 32ð1ð1αα1Þ1=2 1 DSm;1 2 h  i : N DSm1 ¼ ¼ C2  2=5 ρf DSm;2 0:0447C ð1  αÞN vr þ 0:0258N α2

We

6. Evaluation of interfacial area correlation 6.1. Determination of constants The proposed model has been benchmarked against the data sets of Schlegel et al. (2012, 2014) in order to determine the values of the various constants. In order to do so a least-squares regression was used. The error used in the regression was the Root-Mean-Square (RMS) relative error for the Group 2 Sauter mean diameter, Group 1 Sauter mean diameter, and total interfacial area concentration, calculated as

ε¼

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u N    N  N  P P u P DSm2;pred;i  DSm2;meas;i 2 DSm1;pred;i  DSm1;meas;i 2 ai;pred;i  ai;meas;i 2 u þ þ DSm2;meas;i DSm1;meas;i ai;meas;i ti ¼ 1 i¼1 i¼1 3N

:

ð41Þ

Predicted Bubble Size, DSm1,Pred [m]

The values of the constants determined using this method were C1 ¼0.710, C2 ¼ 15.1, and C3 ¼4.30. A comparison of the bench-

0.006 +25%

0.004

0.002

-25%

D = 0.152 m j f = 0.25 m/s j f = 0.50 m/s j f = 1.00 m/s j f = 2.00 m/s

ð220 Þ

ρg

0.000 0.000 0.002 0.004 0.006 Measured Bubble Size, DSm1,Meas [m]

2 bubbles (Shen et al., 2014). The current value of C3 represents an average for both types of injection. Improvements in the correlation prediction could be realized by provided one value of C3 for bubbly flow injection, and a different value for cap-bubbly flow injection. Unfortunately such an approach is not warranted at this time due to the limited amount of experimental data available. In addition to these concerns, many simplifying assumptions were included in the derivation. Given these concerns, this amount of error still represents an improvement over previous models. The variation of the Group 1 bubble size with void fraction is shown in Fig. 10 for various liquid velocities. The plot also shows the limit on the Group 1 bubble size as well as the Group 1 Sauter mean diameters predicted by RELAP and TRACE. The plot reveals some limitations for the new model. This model was developed for cases where the Group 2 void fraction is significant. Thus there is a physically unrealistic trend for the Group 1 Sauter mean diameter as the void fraction becomes small, in the bubbly flow regime. This trend is due to the assumption that Group 1 bubbles are produced by

Predicted Bubble Size, DSm2,Pred [m]

1

0.06

D = 0.152 m jf = 0.25 m/s

+25%

jf = 0.50 m/s

0.04

jf = 1.00 m/s jf = 2.00 m/s

-25%

0.02

0.00 0.00 0.02 0.04 0.06 Measured Bubble Size, DSm2,Meas [m]

Fig. 9. Comparison of predicted and measured Sauter mean diameters.

182

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

10

jf = 0.25 m/s

RELAP

jf = 0.5 m/s

TRACE

Sauter Mean Diameter, Dsm,2[m]

Sauter Mean Diameter, Dsm,1[m]

0.02

jf = 1.0 m/s DSm = 4*Lo

0.01

0.00 0.2

D = 0.152 m P = 180 kPa

DSm = 40*Lo

-1

10

DSm = 4*Lo

-2

0.3

0.4 0.5 0.6 0.7 Void Fraction, α [-]

0.8

10

0.9

0.2

0.3

0.4 0.5 0.6 0.7 0.8 Void Fraction, α [-]

0.9

Interfacial Area Concentration, ai [m-1]

104 TRACE correlation RELAP correlation

103

102

Model Data (Schlegel et al, 2014)

jf = 0.25 m/s jf = 0.5 m/s jf = 1.0 m/s jf = 2.0 m/s D = 0.152 m, P = 180 kPa

101 0.0

0.2

0.4 0.6 0.8 Void Fraction, α [-]

1.0

Ratio of Interfacial Area Concentration Prediction to Measurement, ai,pred / ai, meas [-]

Fig. 10. Bubble size prediction as a function of void fraction.

10 1

Schlegel, et al (2014) D = 0.152 m D = 0.203 m D = 0.304 m

Schlegel, et al (2012) D = 0.152 m D = 0.203 m

+25%

10

0

-25%

10 -1 0.0

0.2

0.4 0.6 0.8 Void Fraction, α [-]

1.0

Fig. 11. Overall interfacial area concentration prediction accuracy.

shearing-off from large Group 2 bubbles. Therefore in order to accurately predict the interfacial area concentration at low void fractions a bubbly flow model such as that of Hibiki and Ishii (2002) should be used to predict the Group 1 Sauter mean diameter. The new model should be phased in using a transition model, such as a linear ramp or exponential transition, as the Group 2 void fraction increases from 0 to 0.1. Over the applicable range the Sauter mean diameter predicted by the new model for Group 1 is above 2 mm but always remains below the upper limit on Group 1 bubble size. Fig. 10 also shows a comparison of the predicted Group 2 bubble size as a function of void fraction for the new model as well as the RELAP and TRACE approaches. In general Group 2 bubbles do not appear at very low void fractions, however under some special conditions they may occur. Thus the model behavior is not plotted for void fractions below 0.2. The shape of the curve for intermediate void fractions reflects the competition between breakup and coalescence effects and remains in the region between the minimum and maximum cap bubble sizes. Above void fractions of 0.9 the flow transitions to annular flow and the model is no longer valid, so again the model is not plotted in this region. An evaluation of the prediction of the total interfacial area concentration can be found in Fig. 11(a) and (b). In Fig. 11(a) the total interfacial area concentration is plotted against void fraction with the different liquid velocity conditions noted. The RELAP and TRACE correlations plotted are for liquid velocity of 0.5 m/s. The

figure shows that the proposed correlation is able to account for the differences between the various liquid flow conditions and is able to accurately capture the change in the interfacial area concentration as the gas flow rate, and therefore void fraction, is increased. This indicates a very strong improvement over the previous correlations used in RELAP and TRACE. A more general comparison, including the data from all of the pipe sizes used by Schlegel et al. (2014) as well as the data of Schlegel et al. (2012), is shown in Fig. 11(b). As in the previous figures, there is a relatively large random error. Most of the data is predicted to within 25%, with an RMS error of 22.6% and a bias of  4.29%. This is consistent with the slightly positive bias in the Sauter mean diameter models, and as stated previously is a reasonable accuracy given the complexity of the phenomenon and the experimental uncertainty. These figures indicate that the new model is capable of predicting interfacial area concentration in large diameter pipes, capturing the effect of various changes in the flow conditions on the hydrodynamics of the flow and improving the ability of existing models to accurately evaluate two-phase flows. The current model cannot be used for small diameter channels. That would require additional assumptions regarding the shape of the interface and the hydrodynamics of the flow. As the data used to benchmark the model was limited to a very narrow pressure range, it has not been confirmed that the results are applicable to high-pressure flows. Additional analysis based on

1000

D = 0.152 m, jf = 0.50 m/s Data Model jg = 0.28 m/s

800

jg = 0.50 m/s jg = 0.95 m/s

600

jg = 4.05 m/s jg = 5.77 m/s jg = 8.55 m/s

400 200 0

0

Interfacial Area Concentration, ai [1/m]

Fig. 12. Axial flow development for 0.152 m diameter channel, jf ¼0.50 m/s.

jg = 0.24 m/s

Interfacial Area Concentration, ai [1/m]

D = 0.304 m, jf = 0.35 m/s Data Model jg = 0.70 m/s jg = 1.14 m/s jg = 1.57 m/s

400

D = 0.304 m, jf = 0.50 m/s Data Model jg = 0.50 m/s

200

jg = 0.84 m/s jg = 1.28 m/s jg = 2.14 m/s

0

5 10 Dimensionless Length, z/D [-]

15

Fig. 16. Axial flow development for 0.304 m diameter channel.

6.2. Applicability for developing flow

jg = 0.52 m/s

800

jg = 0.88 m/s jg = 1.59 m/s

600

jg = 2.95 m/s jg = 4.24 m/s jg = 9.56 m/s

400 200 0

5 10 15 20 25 30 Dimensionless Length, z/D [-]

Fig. 13. Axial flow development for 0.152 m diameter channel, jf ¼1.0 m/s.

1600

D = 0.152 m, jf = 2.00 m/s Data Model jg = 0.45 m/s

1200

jg = 0.94 m/s jg = 1.48 m/s jg = 2.69 m/s jg = 4.35 m/s

800

jg = 5.48 m/s jg = 9.14 m/s

400

0

0

5

10

15

20

25

30

Dimensionless Length, z/D [-] Fig. 14. Axial flow development for 0.152 m diameter channel, jf ¼2.0 m/s.

Interfacial Area Concentration, ai [1/m]

183

D = 0.152 m, jf = 1.00 m/s Data Model

1000

0

600

0

5 10 15 20 25 30 Dimensionless Length, z/D [-]

1200

Interfacial Area Concentration, ai [1/m]

Interfacial Area Concentration, ai [1/m]

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

1000

D = 0.203 m, j = 0.50 m/s Data Model

f

j = 0.68 m/s

800

g

j = 2.20 m/s g

j = 3.56 m/s g

600

D = 0.203 m, j = 1.00 m/s f

Data Model

400

j = 0.37 m/s g

j = 1.36 m/s g

j = 2.95 m/s

200

g

j = 2.15 m/s g

0

0

5

10

15

20

25

Dimensionless Length, z/D [-] Fig. 15. Axial flow development for 0.203 m diameter channel.

experimental data at high pressures, and possibly a revision to the constants determined from the benchmark, is necessary before this model can be considered for a wide range of pressure conditions.

The ability of the correlation to predict the development of two-phase flows along the flow direction is shown in Figs. 12–16. Along the flow direction, the gas volumetric flux continually increases due to expansion of the gas phase as the pressure drops. It should be noted that the correlation assumes that the gas exchange between the two bubble groups is in equilibrium. That is, it assumes that the flow is fully developed. The data used to benchmark the model was restricted to the data at the upper measurement locations, while the inlet measurements were not used since the interfacial area formulation assumed fully-developed flow. The ability of the correlation to predict the inlet conditions and flow development should therefore be confirmed. It was not expected that the correlation would predict the flow development near the entrance to the test section well, as flow at the inlet to the test section is not in equilibrium and is therefore undergoing rapid changes. After equilibrium is reached, the correlation should be capable of predicting the trends in the flow development. The axial change in interfacial area concentration is predicted as follows. First, the inlet pressure and the superficial velocities for each phase are provided to a computer script. Using that data, the script calculates the void fraction, interfacial area concentration, phase velocities, and so on. Then the script moves on the next spatial step in explicit fashion, using the calculated values at the first spatial step to calculate the gravitational and friction pressure drop in a first-order upwind numerical scheme. The friction pressure drop is calculated using the Lockhart and Martinelli (1949) model. The results of this calculation are compared with the experimental data in Figs. 12–16. The data collected for liquid superficial velocity of 0.25 m/s is not shown in the figures. At such a low liquid velocity, the void fraction calculations in this region may contain significant errors which will then carry through to the interfacial area concentration calculations (Shen et al., 2014). Thus these figures only show the data for higher liquid velocity conditions. As expected and explained above, the predictions do improve significantly as the liquid superficial velocity increases. This is mostly attributable to improvements in the prediction accuracy of the void fraction using the drift-flux model, as the distribution parameter becomes more significant than the drift velocity at higher total volumetric fluxes. At liquid velocities higher than 1 m/ s the results are very good, while at liquid velocities of 0.5 m/s and below the results are very scattered. The figures also show that for very high gas velocity conditions the correlation predictions, while still within 20%, are less accurate than predictions for lower gas velocity conditions. Under these conditions the residence time of the fluid mixture in the test section is as low as 0.3 s. This is not enough time for the mixture to

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

Interfacial Area Concentration, ai [m-1]

104 TRACE correlation RELAP correlation

103

102 Model Data (Yang et al, 2012) j = 1.0 m/s f

j = 0.5 m/s f

j = 0.25 m/s f

101 0.0

0.2

0.4 0.6 0.8 Void Fraction, α [-]

1.0

Ratio of Interfacial Area Concentration Prediction to Measurement, ai,pred / ai, meas [-]

184

101 Yang, et al. (2012)

+25%

10 0 -25%

-1

10

0.0

0.2

0.4 0.6 0.8 Void Fraction, α [-]

1.0

Fig. 17. Comparison of model with rod bundle data.

reach equilibrium between the two bubble groups. Because of this some degradation is expected, but compared to the experimental error (10%) and the overall bias and error of the correlation (  4.29% and 22.6%, respectively) this degradation is minor. The overall conclusion based on these figures is that the correlation can predict changes in the interfacial area concentration due to gas expansion, provided that the correlation is not used too close to the flow channel inlet and that sufficiently accurate void fraction models are used. However like many other algebraic correlations it was developed based on the assumption of equilibrium conditions. Therefore care must be used when applying the correlation to developing flows with small values of z/D, such that the flow is not able to reach equilibrium. 6.3. Applicability for rod bundles Flows in rod bundles are a major concern in nuclear reactor systems and various heat exchanger designs. There are a number of similarities between flows in rod bundles and flows in large diameter channels. In rod bundles, large Group 2 bubbles might occupy an entire subchannel but very rarely occupy the entire flow area. This is especially true for Pressurized Water Reactor (PWR) rod bundles, where two-phase flows occur during most nuclear power plant accident scenarios. In this sense, the physical phenomena governing the coalescence and disintegration of bubbles in the two geometries may be very similar. In order to test the performance of the proposed correlation for rod bundles, the correlation results have been compared to data collected in an 8  8 simulated BWR rod bundle geometry by Yang et al., 2013. The results of the comparison are shown in Fig. 17. Fig. 17(b) shows that the model is able to predict most of the data at lower void fractions quite well, with an average RMS error of 23.6% and a bias of 12.0%. Fig. 17(a) shows that the data for the rod bundle is much closer than the data for large diameter channels to the predictions given by the models used in TRACE, but that the new correlation results in improved predictions for the interfacial area concentration. The performance of the correlation begins to degrade as the flow transitions out of cap-bubbly flow and into churn-turbulent flow. This occurs at void fractions of 0.4 to 0.5. This degradation occurs because in this range, cap bubbles begin to occupy multiple subchannels. This puts additional constraints on the shape and size of Group 2 bubbles that are not accounted for in the development of the correlation. The results may be improved by modifying the coefficient in the Group

2 Sauter mean diameter model for this specific geometry, however, the existing database is insufficient to justify such a change. The results might also be improved by implementing a turbulence model that accounts for the effects of the spacer grid, as the measurement in the region following the spacer grids is often the measurement with the highest prediction error.

7. Conclusions In this paper, a new algebraic correlation has been presented for the prediction of interfacial area concentration in two-phase flows based on a two-bubble-group method. Some specific conclusions of this effort include:


Constants were benchmarked based on the data of Schlegel















et al. (2014). The resulting models have an average bias and prediction error of  4.29% and 22.6% for the total interfacial area concentration, in comparison with the data collected by Schlegel et al. (2012, 2014). The correlation is validated for air–water flows at pressures of 100 kPa to 300 kPa, liquid velocities of 0.25 to 2 m/s, and void fractions of 0.2 to 0.95. The correlation is validated for pipes with diameters from 0.152 m to 0.304 m, but is expected to perform well even in pipes with much larger diameters. Performance could be improved by providing separate values of C3 for bubbly flow and cap bubbly flow injection, if sufficient data was available for validation. The model also performs reasonably well for rod bundle geometries at void fractions lower than 0.5. Performance could be improved by modifying the coefficient C3, if sufficient data was available for validation. Comparison with the models included for high void fraction flows in TRACE and RELAP show that the new model greatly improves the ability to predict interfacial area concentration in cap-bubbly and churn turbulent flows, and results in minor improvements for flows with void fractions between 0.1 and 0.3. The correlation is able to capture various physical flow behaviors, including the transition to annular flow and the effect of flow conditions on the bubble sizes. In addition to calculating the total interfacial area concentration this model allows the calculation of the void fractions and Sauter mean diameters of each bubble group, so that this model could be fit into a two-group two-fluid model calculation.

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186


The collection of additional databases including area-averaged interfacial area concentration would allow improved benchmarking and expansion of the range of conditions that can be predicted using this correlation, as well as the accuracy of the correlation.

TI WE

185

turbulent impact Wake entrainment

References NomenclatureLatin characters ai C C0 CD Dh DSm f tp G g j Lo N μf N DSm N Rem N Reε N vr N We N Wec v vgj

interfacial area concentration [m  1] model constant [-] distribution parameter [-] drag coefficient [-] hydraulic diameter [m] Sauter mean diameter [m] two-phase friction factor [-] mass flux [kg  m  2 s  1] gravitational acceleration [m s  2] volumetric flux [m s  1] Laplace length scale [m] viscosity number [-] Sauter mean diameter ratio [-] mixture Reynolds number [-] dissipation Reynolds number [-] relative velocity ratio [-] Weber number [-] critical Weber number [-] velocity [m s  1] drift velocity [m s  1]

Greek characters

α Δρ ε εb εBI εsh φ ν ρ σ

void fraction [-] density difference [kg m  3] RMS error [-] error bias [-] bubble-induced dissipation [m2 s  3] shear induced dissipation [m2 s  3] interfacial area source/sink [m  1 s  1] kinematic viscosity [m2 s  1] density [kg m  3] surface tension [N m  1]

Superscript/subscript þ,n 1 2 1; 2 2; 1 base bs c crit f g gs L m max r RC sa SO

dimensionless value value for group 1 value for group 2 interaction between bubble groups interaction between bubble groups base value value at the bubbly-slug transition critical value critical value value for the liquid value for the gas value in the liquid slug limiting value mixture value maximum value relative value random collision value at the slug-annular transition shearing-off

Akita, K., Yoshida, F., 1974. Bubble size, interfacial area, and liquid-phase mass transfer coefficient in bubble columns. Ind. Eng. Chem. Process Des. Develop 13, 84–91. Brooks, C.S., Paranjape, S., Ozar, B., Hibiki, T., Ishii, M., 2012. Two-group drift-flux model for closure of the modified two-fluid model. Int. J. Heat Fluid Flow 37, 196–208. Fu, X.Y., Ishii, M., 2003a. Two-group interfacial area transport in vertical air–water flow: I. Mechanistic model. Nucl. Eng. Des. 219, 143–168. Fu, X.Y., Ishii, M., 2003b. Two-group interfacial area transport in vertical air–water flow: II. Model evaluation. Nucl. Eng. Des. 219, 169–190. Hazuku, T., Takamasa, T., Hibiki, T., Ishii, M., 2007. Interfacial area concentration in annular flow. Int. J. Heat Mass Transfer 50, 2986–2995. Hibiki, T., Ishii, M., 2000a. Two-group interfacial area transport equations at bubblyto-slug flow transition. Nucl. Eng. Des. 202, 39–76. Hibiki, T., Ishii, M., 2000b. One-group interfacial area transport of bubbly flows in vertical round tubes. Int. J. Heat Mass Transfer 43, 2711–2726. Hibiki, T., Ishii, M., 2001a. Interfacial area concentration in steady fully-developed bubbly flow. Int. J. Heat Mass Transfer 44, 3443–3461. Hibiki, T., Ishii, M., 2001b. Effect of inlet geometry on hot-leg U-bend two-phase natural circulation in a loop with a large diameter pipe. Nucl. Eng. Des. 203, 209–228. Hibiki, T., Ishii, M., 2002. Interfacial area concentration of bubbly flow systems. Chem. Eng. Sci. 57, 3967–3977. Hibiki, T., Ishii, M., 2003. One-dimensional drift-flux model for two-phase flow in a large diameter pipe. Int. J. Heat Mass Transfer 46, 1773–1790. Hibiki, T., Lee, T.H., Lee, J.Y., Ishii, M., 2006. Interfacial area concentration in boiling bubbly flow systems. Chem. Eng. Sci. 61, 7979–7990. Huh, B.G., Euh, D.J., Yoon, H.Y., Yun, B.J., Song, C.H., Chung, C.H., 2006. Mechanistic study for the interfacial area transport phenomena in an air–water flow condition by using fine-size bubble group model. Int. J. Heat and Mass Transfer 49, 4033–4042. Ishii, M., Hibiki, T., 2010. Thermo-Fluid Dynamics of Two-Phase Flow. Springer, New York. Ishii, M., 1977. One dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow regimes. ANL, 77 47. Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly, droplet or particulate flows. AIChE J. 25, 843–855. Ishii, M., Mishima, K., 1984. Two-fluid model and hydrodynamic constitutive relations. Nucl. Eng. Des. 82, 107–126. Kataoka, I., Ishii, M., 1987. Drift-flux model for large diameter pipe and new correlation for pool void fraction. Int. J. Heat Mass Transfer 30, 1927–1939. Lehr, F., Mewes, D., 2001. A transport equation for the interfacial area density applied to bubble columns. Chem. Eng. Sci. 56, 1159–1166. Liao, Y., Lucas, D., Krepper, E., Schmidtke, M., 2011. Development of a generalized coalescence and breakup closure for the inhomogeneous MUSIG model. Nucl. Eng. Des. 241, 1024–1033. Lin, C., Hibiki, T., 2014. Databases of interfacial area concentration in gas–liquid two-phase flow. Prog. Nucl. Energy 74, 91–102. Lockhart, R.W., Martinelli, R.C., 1949. Proposed correlation of data for isothermal two-phase, two-component flow in pipes. Chem. Eng. Prog. 45, 39–48. Lucas, D., Beyer, M., Szalinski, L., Schutz, P., 2010. A new database on the evolution of air–water flows along a large vertical pipe. Int. J. Therm. Sci. 49, 664–674. Ohnuki, A., Akimoto, H., 2000. Experimental study on transition of flow pattern and phase distribution in upward air–water two-phase flow along a large vertical pipe. Int. J. Multiphase Flow 26, 367–386. Ohnuki, A., Akimoto, H., 2001. Model development for bubble turbulent diffusion and bubble diameter in large vertical pipes. J. Nucl. Sci. Technol. 38, 1074–1080. Ozar, B., Jeong, J.J., Dixit, A., Julia, J.E., Hibiki, T., Ishii, M., 2008. Flow structure of gas–liquid two-phase flow in an annulus. Chem. Eng. Sci. 63, 3998–4011. Ozar, B., Dixit, A., Chen, S.W., Hibiki, T., Ishii, M., 2012. Interfacial area concentration in gas–liquid bubbly to churn-turbulent flow regime. Int. J. Heat Fluid Flow 38, 168–179. Prasser, H.M., Beyer, M., Carl, H., Gregor, S., Lucas, D., Pietruske, H., Schutz, P., Weiss, F.P., 2007. Evolution of the structure of a gas–liquid two-phase flow in a large diameter pipe. Nucl. Eng. Des. 237, 1848–1861. Sari, S., Ergun, S., Barik, M., Kocar, C., Sokmen, C.N., 2009. Modeling of isothermal bubbly flow with interfacial area transport equation and bubble number density approach. Ann. Nuclear Energy 36, 222–232. Schlegel, J.P., Hibiki, T., Ishii, M., 2010. Development of a comprehensive set of driftflux constitutive models for pipes of various hydraulic diameters. Prog. Nucl. Energy 52, 666–677. Schlegel, J.P., Miwa, S., Chen, S., Hibiki, T., Ishii, M., 2012. Experimental study of twophase flow structure in large diameter pipes. Exp. Therm. Fluid Sci. 41, 12–22. Schlegel, J.P., Macke, C., Hibiki, T., Ishii, M., 2013. Modified distribution parameter for churn-turbulent flow in large diameter channels. Nucl. Eng. Des. 263, 138–150.

186

J.P. Schlegel, T. Hibiki / Chemical Engineering Science 131 (2015) 172–186

Schlegel, J.P., Sharma, S., Cuenca, R.M., Hibiki, T., Ishii, M., 2014. Local flow structure beyond bubbly flow in large diameter channels. Int. J. Heat Fluid Flow 47, 42–56. Serizawa, A., Kataoka, I., 1990. Turbulence suppression in bubbly two-phase flow. Nucl. Eng. Des. 122, 1–16. Shen, X., Schlegel, J.P., Chen, S.W., Rassame, S., Griffiths, M.J., Hibiki, T., Ishii, M., 2014. Prediction of Void Fraction in Large Diameter Pipes Using Drift-Flux Models. Springer Briefs on Frontiers and Progress in Multiphase Flow I, Chapter 2. Springer, New York. Smith, T.R., Schlegel, J.P., Hibiki, T., Ishii, M., 2012a. Two-phase flow structure in large diameter pipes. Int. J. Heat Fluid Flow 33, 156–167. Smith, T.R., Schlegel, J.P., Hibiki, T., Ishii, M., 2012b. Mechanistic modeling of interfacial area transport in large diameter pipes. Int. J. Multiphase Flow 47, 1–16.

Spore, J.W., Jolly-Woodruff, S.J., Knight, T.K., Lin, J.C., Nelson, R.A., Pasamehmetoglu, K.O., Steinke, R.G., Unal, C., 1993. TRAC-PF1/MOD2, vol. 1: Theory Manual. LA-12031-M, vol. I, NUREG/CR-5673, USA. Sun, X., et al., 2003. Modified two-fluid model for the two-group interfacial area transport equation. Ann. Nucl. Energy 30, 1601–1622. Thermal Hydraulics Group, 1998. Relap5/Mod3 Code Manual. Vol. I: Code Structure, System Models and Solution Methods. RELAP5/MOD3.2.2Beta, USA. Yang, X., Schlegel, J.P., Liu, Y., Paranjape, S., Hibiki, T., Ishii, M., 2013. Experimental study of interfacial area transport in air–water two-phase flow in a scaled 8  8 BWR rod bundle. Int. J. Multiphase Flow 50, 16–32. Yao, W., Morel, C., 2004. Volumetric interfacial area prediction in upward bubbly two-phase flow. Int. J. Heat Mass Transfer 47, 307–328.