Abstract
The one-dimensional two-fluid model approach has been traditionally used in thermal-hydraulics codes for the analysis of transients and accidents in water–cooled nuclear power plants. This paper investigates the performance of RELAP5/MOD3 predicting vertical upward bubbly flow at low velocity conditions. For bubbly flow and vertical pipes, this code applies the drift-velocity approach, showing important discrepancies with the experiments compared. Then, we use a classical formulation of the drag coefficient approach to evaluate the performance of both approaches. This is based on the critical Weber criteria and includes several assumptions for the calculation of the interfacial area and bubble size that are evaluated in this work. A more accurate drag coefficient approach is proposed and implemented in RELAP5/MOD3. Instead of using the Weber criteria, the bubble size distribution is directly considered. This allows the calculation of the interfacial area directly from the definition of Sauter mean diameter of a distribution. The results show that only the proposed approach was able to predict all the flow characteristics, in particular the bubble size and interfacial area concentration. Finally, the computational results are analyzed and validated with cross-section area average measurements of void fraction, dispersed phase velocity, bubble size, and interfacial area concentration.
1. Introduction
Two-phase flow phenomena have been an object of study during several decades with a great impact in nuclear field. From the reactor to the turbines, one can find a wide variety of systems where two-phase flow plays a main role: BWR core, secondary loop, or reactor heat removal system (RHRS) are examples of two-phase flow components. In such cases, two-phase flow is present in normal operating conditions, but also in specific situations, like instabilities events, loss-of-coolant accidents, or refueling. The previous cases imply different conditions of pressure, temperature, or mass flow.
This broad range of situations is considered in one-dimensional thermal-hydraulics codes to set the appropriate flow regime in each situation. They are based on the two-fluid model [1], where averaged Navier-Stokes equations are solved for each phase including momentum, energy, and continuity equations. Then, one can account for the interaction terms between phases, to consider the mass transfer, momentum, and energy at the interface. The interfacial momentum term differs depending on which flow regime is working. The proper regime is selected according to a flow regime map and the velocities of each phase. Different flow regime maps have been proposed by different authors [2, 3]. This paper investigates the performance of RELAP5/MOD3 predicting the results of experiments in a vertical upward bubbly flow for low velocity conditions. Bubbly flow at those conditions can be found in pressurizers, reactor pools, or refueling operations. To investigate in depth the bubbly flow behavior and the one-dimensional modeling, we perform experiments in a vertical pipe of around 6 meters of length along which an air-water fluid in bubbly regime moves upwards in adiabatic conditions and atmospheric pressure.
In the one-dimensional two-fluid model (1D TFM) the drag term plays the main role in the interfacial momentum transfer. System codes use different approaches to model the interfacial drag force depending on the flow regime. Two approaches are usually used to define the interfacial drag force: the drift-velocity approach (DVA) and the drag coefficient approach (DCA).
RELAP5/MOD3 uses DVA for bubbly flow in vertical pipes and DCA was used in the previous version RELAP5/MOD2. The drift models, although simpler than the drag coefficient approach, are usually only valid in the range of applicability for which they were obtained as they depend on flow and geometry.
In this work we make use of DCA in RELAP5/MOD3, by modifying the code. The drag force calculated with DCA relies on correlations that are defined traditionally as a function of Reynolds and/or Eötvös numbers, making the calculation of the drag term more general than with the drift-velocity.
However, the use of DCA incorporates a set of assumptions to calculate the drag term. We propose a more rigorous version of the drag coefficient approach () that has been implemented to evaluate the influence of these assumptions, consisting of the following:(i)A drag coefficient correlation that takes into account the effect of the bubble shape through the Eötvös number and the effect of the contaminants present in the system used.(ii)Bubble size distribution (BSD) consideration by means of its statistical parameters including the axial evolution due to the gas expansion.(iii)Interfacial area calculated directly from the definition of the Sauter mean diameter of the BSD.
In summary, three drag coefficient approaches are used: a drift-velocity approach named DVA, an existing drag coefficient approach DCA, and the proposed drag coefficient approach (see Figure 1).
These approaches are assessed comparing the influence of the simplifications and assumptions on the results. Then, a validation of the cross-section average void fraction, dispersed phase velocity, bubble size, and interfacial area concentration is performed.
This paper is structured as follows. First, in Section 2 a description of the experimental data and the measurements of the bubble size along the pipe are presented. In Section 3, the mathematical formulation of the different approaches is introduced. Section 4 describes the model and setup of the simulation. Section 5 includes a model assessment with different drag approaches. Section 6 shows the validation between computational results and experimental data. Finally, some conclusions are drawn in Section 7.
2. Experimental Facility and Measurements of Bubble Size Distribution
2.1. Description of the Experimental Facility and Techniques
The experimental data used in this paper are based on the facility described in Monrós-Andreu et al. [6, 7] in a pipe of diameter 52 mm and length 5500 mm with a sparger to inject the air flow. Three cases with liquid velocity 0.5 m/s and different gas superficial velocities are investigated: 0.02 m/s, 0.03 m/s, and 0.04 m/s, labeled PW05002, PW05003, and PW05004, respectively. The details about the experimental facility and measurements are described below. The experimental facility is located at the Laboratory of Hydraulics of the Universitat Jaume I and consists of an upward flow experimental loop (Figure 2).
Osmotic water (200–300 μS m−1) was circulated by a centrifugal pump and stored in a 500 L reservoir tank at constant temperature (20°C) by a heat exchanger. The water flow rate introduced in the system was measured by an electromagnetic flow meter (M1000, Badger Meter Inc.). An air flow-meter controller (EL-FLOW 250 lNpm, Bronckhorst Hi-Tech) was used to adjust and measure the gas flow rate. Three axial locations were used for the measurements: (bottom), (middle), and (top). Conductivity probes and LDA measurements were performed in these three ports, while high-speed cameras are located at bottom and top ones.
The measurement system for the dispersed phase consisted of three mounted four-sensor conductivity probes, mechanical traversers, a measurement circuit, a digital high-speed acquisition board, and the software used for signal processing. The four-sensor conductivity probe was attached to the mechanical traverser mounted on a specially designed flange, and it moved along the radial direction of the test section using controlled step motors. The measurement circuit was used to measure the voltage difference between the exposed tip and the grounded terminal. A high-speed acquisition board (National Instrument Crop., SCXI-1325) and a PC were used to acquire the signals of the four-sensor probe, with a control program developed under the LabView (National Instrument Crop.) software environment.
Pressure was measured by pressure transducers with a range 0-1 bar for the two lower ports, and 0–200 mBar for the mid and top ports (all transducers with an 0.1% relative error).
2.2. Bubble Size Distribution Measurements
In order to evaluate the drag coefficient approaches, the BSD was first investigated experimentally at different heights for the different flow conditions by analyzing the images obtained by the high-speed cameras. An example of these image is shown in Figure 3.
Approximately 500 bubbles were manually measured for every port similarly as in Peña-Monferrer et al. [8]. For the measurements, several points in the bubble borders are selected and an ellipse is fitted to the selected points, using a least-squares algorithm that provides both axis and orientation angle. The semiaxes are used to obtain for each bubble. The BSD of the equivalent diameters at = 22.4 and = 98.7 are shown in Figure 4. In the figure, the histogram at both heights and the normal distribution fitted are shown. At first glance, an increase in bubble size between the bottom and the top measurement ports is observed for the three cases as a result of the bubble decompression. The BSD was fitted in the literature for bubbly flow to normal [9], log-normal [10–14], or gamma [15, 16] distributions. For the experiments in this paper, where relatively low superficial gas velocities are given, the bubble size data fitted well to a normal distribution. All the samples passed the Kolmogorov-Smirnov test at 5% significance level. It was applied for PW05002, PW05003, and PW05004. The values for these cases at are 0.26, 0.60, and 0.94 respectively.
3. Mathematical Formulation and Implementation
3.1. Drift-Velocity Approach (DVA)
The relative motion between the phases can be considered through a drift-flux model [17, 18]. In a drift-flux model, the mixture of the phases is solved as a whole. In RELAP5/MOD3, a drift-velocity approach is incorporated into the TFM to describe the interfacial drag force term. The area-averaged interphase drag term is given aswhere is the drag coefficient and the relative velocity between both phases. The drag coefficient is obtained from a balance of the forces in the direction of the flow. It considers the interfacial drag, buoyancy, and pressure drop and applies the assumptions that both phases have equivalent pressure and the action-reaction principle for the interfacial momentum terms [19]. The drag coefficient is given bywhere , , and are the void fraction, gravitational constant, and density. In this equation, the relative velocity is replaced with the relation between local relative velocity and void weighted phase velocities assuming uniform relative velocity:
From the definitions of (2) and (3), the drag coefficient in terms of drift-flux is finally given as
The term refers to the void fraction weighted average drift velocity that depends on the flow geometry. For vertical pipe flows and conditions studied in this work, RELAP5/MOD3 use the Chexal-Lellouche correlation [20, 21]:where is the surface tension. The equation depends on many constants as , , , and among others. This is a generalized correlation that was compared using steam-water, air-water, and refrigerant data on multiple flow configurations. They are ranging from different orientations as vertical, horizontal, or inclined, different geometries as pipes, channels, rod bundle, and flow configurations as cocurrent or countercurrent flow. This correlation, although general, lacks the model specificity that is required for an accurate prediction [22]. Moreover, this approach is not consistent with the TFM and its application is contrary to the field equations solved as noted by Brooks et al. [19].
3.2. Drag Coefficient Approach (DCA)
The drag coefficient approach is based on the general drag interfacial term. This is defined aswhere , , and are the drag coefficient, interfacial area concentration, and difference between area-averaged mean velocities. For RELAP5/MOD3 and bubbly flow, the drag coefficient is based on Ishii and Chawla [5]:
The drag coefficient depends on the flow parameters, and the bubble size. The maximum bubble diameter is calculated from a critical Weber number:where and are the surface tension and bubble diameter, respectively.
RELAP5/MOD3 specifies a value of 10 for bubbles for the [23, 24]. In this equation, is not calculated as the difference between the phase velocities but refers to the velocity difference that gives the maximum bubble size [23]. This velocity is also used to calculate the Reynolds number. The following equation is applied:where is set to 0.005 m for bubbly flow and is the hydraulic diameter.
The bubble diameter is calculated from the maximum bubble diameter with the following assumption:
The interfacial area concentration is then given in terms of the mean bubble diameter [19, 23]:where is the Sauter mean diameter of the distribution related to the bubble diameter assuming a Nukiyama-Tanasawa distribution [23], a distribution for droplet diameter for a spray.
3.3. Drag Coefficient Approach with Specific Drag Closure and Bubble Size Distribution ()
The previous drag coefficient approach contains a set of assumptions that may affect the prediction of the two-phase flow characteristics. In this work we propose a new approach to validate bubbly flow scenarios. The model consists of a drag coefficient that takes into account the bubble size effects. On one hand, it considers the bubble dynamics as a function of the bubble size and shape. On the other hand, the BSD is specified as an inlet boundary condition and is used to directly compute the interfacial area. The evolution of the BSD is calculated to account for the decompression effect on the size.
The drag correlation of Tomiyama et al. [4] for contaminated systems is used and implemented:
The Eötvös number is defined as follows:where is the gravity constant.
This expression includes a region dominated by the Eötvös number. Thus, an appropriate calculation of the drag coefficient requires an accurate representation of the bubble size in the system. The terminal velocity of the bubbles as function of the equivalent diameter using the drag force coefficients of (7) and (12) is compared in Figure 5.
The bubble size at the inlet boundary condition is defined from the measurements of the BSD in the experiments. The differences in the bubble size between two heights in a pipe, excluding the breakup and coalescence mechanisms, are due to the pressure changes. As the experiments are performed at atmospheric pressure, an increase of around 30% can be noted from the bottom to the top ports in the experiments ( to ). For a rigorous implementation, the axial change on the bubble size must be considered. Note that breakup and coalescence are neglected because of the flow conditions. The observations with the high-speed camera confirmed this fact. For other scenarios a one-dimensional approximation of a population balance equation would be required, but for this work this approximation has been preferred for convenience as a first approach.
This change in size can be described by an expansion factor that is related in this work for convenience to the inlet values. At each node we can calculatewhere and are the void fraction at the given node and the void fraction at the inlet, respectively.
If the bubbles change their size by the factor ; this means a proportional increase of the bubble size and it is equivalent to multiplying a random variable by a constant value. Then, the mean or expected value of the BSD is also multiplied by the constant value and the same is applied to the standard deviation:
Then, the BSD can be estimated as a scaled distribution of the BSD at the different heights or nodes. A normal distribution at a given height would have the following statistical parameters related to the inlet:
A mean bubble diameter of the distribution can be defined from the numeric mean diameter definition:
The Sauter mean diameter of the distribution can be calculated from its general definition knowing that the bubble size follows a normal distribution:
The interfacial area concentration is obtained from the definition of Sauter mean diameter giving:
3.4. RELAP5 Implementation
In order to get simulations with DCA and methods, some modifications were performed in RELAP5. Hereafter, we present a brief description of the implementation developed in the code for each approach. In the case of DCA approach, only small changes were necessary. Then, to implement , we used DCA as a basis, including the Tomiyama drag correlation, the definition of the Sauter mean diameter from the statistical parameters of the distribution, and the interfacial area using this Sauter mean diameter.
RELAP5 includes already the DCA model, but it is used only for horizontal pipes. Therefore, in this case the main target was to create a flag for accessing the subroutine where the drag coefficient approach was implemented. A variable was created to use this subroutine in vertical bubbly flow cases.
4. Modeling and Setup
The simulations are undertaken by modeling a pipe, whose length is equal to the experimental section from = 22.4 (inlet) to = 98.7 (outlet), with 99 uniform axial nodes. A scheme of the model of RELAP5/MOD3 and the nodalization are shown in Figure 6.
The flow conditions used for this work are summarized in Table 1. In this table, the mean and standard deviation parameters of the normal distribution fitting the BSD for the inlet of the simulation (see Figure 2) are shown. For those scenarios, the water velocity was fixed to 0.5 m/s.
Boundary conditions are defined by time-dependent volumes at both inlet and outlet, followed by a time-dependent junction at the inlet and a branch at the outlet. The values shown in Table 1 are defined at the corresponding TMDPVOL component. As a result, the values shown at the first node are obtained from the resolution of the governing equations.
In order to simulate noncondensable gases, one has to activate card 110 in the input. This card allows defining one or more (until eight) gases. In this work, only air has been defined. RELAP5/MOD3 changes from single-phase to two-phase critical flow model when the noncondensable quality is greater than 1 × 10−6. From this moment, the gas phase is treated as a mixture of vapor and noncondensable gas [23].
The simulations performed consisted of a null transient of 100 seconds, where time step was fixed to 1 × 10−3 to guarantee stability, satisfying a Courant number lower than unity.
The simulations are performed with DVA, DCA, and the modified version . In this section we show first a comparison of these models and later a validation with experiments using the proposed model . Cross-section averaged experimental values are obtained from the radial profiles to compare the results of the simulations with the experiments.
5. Model Assessment
The different drag approaches are assessed in this section. Only the PW05003 case is presented for the sake of clarity as the other flow conditions exhibits similar flow characteristics. In the experiments, the included error bars correspond to the standard error, obtained through repeated observations.
Figure 7 shows the axial evolution of void fraction. The effect of the gas decompression is noted in the void fraction as a function of the height. DCA and give similar results while DVA shows lower void fraction values with a smoother axial evolution than the drag coefficient approaches. The discrepancies increase with the height, mainly because takes into account the bubble expansion of the distribution, and in consequence it has an impact on the drag coefficient and the interfacial area.
Figure 8 shows the comparison of the dispersed phase velocity. DVA uses the Chexal-Lellouche and gives a higher dispersed phase velocity and consequently the void fraction values shown before are significantly lower. Slightly different trends are noted with DCA and with decreasing values of the velocity along the pipe.
The bubble diameter is calculated for DCA and (see Figure 9) as they are based on the drag coefficient approach and the size of the bubble is involved. Note that a proper calculation of the bubble diameter could be required to take into account, for example, heat transfer, breakup, or coalescence phenomena. Significant discrepancies in the mean bubble size for both approaches are shown in the figure.
The interfacial area concentration is compared, in turn for both approaches. An accurate prediction of this variable is required to compute the drag term or when heat transfer plays an important role as the simulations in nuclear installations. Figure 10 shows the comparison for these cases.
Note that the values of the interfacial area are relatively close for both approaches. The interfacial area equation for DCA is based on several assumptions but for is calculated directly from the Sauter mean diameter of the BSD. For DCA, the bubble diameter is underestimated but the interfacial area modeling compensates it, obtaining reasonable results for the interfacial area. Using DCA, a mean bubble size similar as the one obtained in the experiments would result in an underestimation of the interfacial area concentration close to 50%. These observations were also appreciated for PW05002 and PW05004 cases. This unrealistic equilibrium between bubble size and interfacial area, for our experiments and using DCA, can be explained by four main factors:(i)The use of the Nukiyama-Tanasawa size distribution.(ii)The criteria to determine the maximum bubble size from a critical Weber number.(iii)The assumption of obtaining the bubble diameter as half of the maximum diameter.(iv)The calculation of the relative velocity between phases with (9).
6. Validation and Discussion
The previous section showed the differences existing when using the different approaches. It demonstrated for a given scenario that void fraction, dispersed phase velocity, mean bubble size, and interfacial area concentration can vary widely if a proper representation of the drag force and bubble size is not considered in the simulation. For instance, common models as DVA or DCA are not able to predict altogether the variables checked for this scenario due to the assumptions included. However, the proposed drag coefficient approach, , was able to predict the flow characteristics evaluated. In this section this model is used to analyze and validate the PW05002, PW05003, and PW05004 scenarios described in Table 1.
The mean bubble size and its axial evolution are shown in Figure 11. Bigger bubble sizes are noted for PW05004 as higher gas flow rates through the sparger could result in an increasing diameter. However, a linear relation is not observed between the mean bubble size and the gas superficial velocity for the three cases. It could be explained due to the mechanism described by Kazakis et al. [13] where as the gas flow rate increases, more sparger pores are activated and hence more bubbles are formed. For higher values, larger bubbles can be produced from the activated pores or eventually if smaller pore sizes exist new smaller bubbles will appear.
The void fraction profiles are compared in Figure 12. The computational results match the experiments accurately along the pipe. The results at the top measurement port are well predicted for the three cases despite its nonlinear evolution. This effect is more pronounced as the gas flow rate increases.
In Figure 13 the validation is done for the dispersed phase velocity. The results are similar to the experiments in both magnitude and trend. The drag coefficients are obtained from experiments for single bubbles and therefore the influence that the bubbles have with each other is not taken into account. Therefore, the dispersed phase velocity of the system could be different with these considerations.
Finally, the interfacial area concentration of the simulation is compared with the experiments in Figure 14. The experimental values are obtained with the definition of interfacial area (19) using the measures of void fraction and BSD.
7. Conclusions
A 1D TFM was used to simulate bubbly flow in adiabatic air-water upward bubbly flow to compare the results with experimental data at low velocity conditions. RELAP5/MOD3 was used to simulate these scenarios. This system code uses the drift-velocity approach (DVA) by default for the conditions tested. The differences with the experiments were considerable using this approach. Then the drag coefficient approach (DCA) used in RELAP5/MOD2 for bubbly flow and vertical pipes was incorporated giving more reasonable results. However, the simulation with DCA was not able to predict all the variables compared, in particular, the mean bubble diameter and the relation with the interfacial area.
A modified version of the drag coefficient approach implemented was proposed (). This included a proper drag force for this scenario, the size effects in the drag force, bubble size distribution (BSD) with axial evolution, and direct calculation of the interfacial area from the definition of the Sauter mean diameter. As a result of this implementation more accurate results in terms of magnitude and trend are obtained overall with regard to DCA. The validation performed with this model showed a good agreement with the experiment for several variables as axial evolution of void fraction, dispersed phase velocity, mean bubble size, and interfacial area. From this study the following conclusions are drawn:(i)The Chexal-Lellouche drift correlation fails in predicting bubbly flow in vertical pipes at low liquid velocities.(ii)DCA predicts relatively well the interfacial area at expenses of underestimating the mean bubble size close to 50%.(iii)The proposed drag coefficient approach is able to reproduce all the variables as the size distribution is considered and interfacial area is calculated directly from the size distribution.
This demonstrates that considering the BSD and its effects is required as shown with the simulations with . It is crucial to obtain a generic model able to predict the interfacial area for different scenarios. In addition, a mean size distribution in accordance with the interfacial area can be calculated. This can be important for the development of submodels for heat transfer, regime transition, or, for example, chemical models when applying 1D TFM to other industries. While the present study was focused on investigating the different drag approaches using a system code, future investigations will incorporate the study of high velocity conditions similar to the present in nuclear reactors where breakup and coalescence take place. These investigations should include 1D TFM modeling using population balance equations.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors sincerely thank the “Plan Nacional de I+D+i” for funding the Projects MODEXFLAT ENE2013-48565-C2-1-P, ENE2013-48565-C2-2-P, and NUC-MULTPHYS ENE2012-34585.