CFD modelling of two-phase discharge from a stratified region through a small side branch

Two-phase discharge from headers has relevance to many industrial applications, particularly in nuclear-safety analysis during postulated loss-of-coolant accidents. The main objective of this study is to assess the ability of commercial CFD codes in predicting the flow phenomena and flow parameters associated with two-phase discharge from a stratified region through a small side branch. Results were obtained for the critical heights at the onsets of gas and liquid entrainment, as well as the mass flow rate and quality during two-phase discharge. All results are in good agreement with existing experimental data. However, the computation time required to obtain these results was found to be excessive.


Introduction
Two-phase discharge from a stratified region through single or multiple branches has relevance to many industrial applications. Examples of these applications include the flow through small breaks in horizontal pipes during postulated loss-of-coolant accidents in nuclear reactors, the flow distribution in multichannel heat exchangers, and the flow distribution in two-phase headers and manifolds.
For the case of discharge from a stratified two-phase region through a single horizontal branch, Zuber [1] pointed out that two distinct flow phenomena may occur depending on the location of the gas-liquid interface relative to the branch. If the branch is located above the horizontal interface, liquid can be entrained in the predominantly gas flow, and this phenomenon is referred to as "liquid entrainment". On the other hand, if the branch is located below the horizontal interface, gas can be pulled through the branch (gas entrainment phenomenon). When the position of the interface relative to the branch is somewhere between the positions for the onset of liquid entrainment (OLE) and the onset of gas entrainment (OGE), a two-phase mixture flows into the branch.
Several experimental (e.g., Micaelli and Memponteil [2];Yonomoto and Tasaka [3]; Hassan et al. [4]) and theoretical (e.g., Craya [5]; Soliman and Sims [6]; Saleh et al. [7]) investigations were reported for the determination of the critical heights of the interface at the OLE and OGE. These studies produced the following formulation for both critical heights: where, h is the critical height for OLE or OGE (measured from the centerline of the branch), d is the branch diameter, Fr is Froude number, and C 1 and C 2 are coefficients whose values depend on the type of the phenomenon. The theoretical studies were in general good agreement with the experimental data.
For the condition where the interface is located between the two critical heights, experimental data were obtained for the mass flow rate and quality of the twophase mixture exiting through the branch (Yonomoto and Tasaka [3]; Hassan et al. [4]). Empirical correlations were obtained from these data (e.g., [4]); however, there have been no theoretical studies, analytical or numerical, for predicting these data. The objective of the present study is to develop numerical models for predicting the critical heights at the OLE and OGE, as well as the mass flow rate and quality in the two-phase region and to validate the numerical results in terms of magnitude and trend against existing experimental data.

Numerical model
The geometry under consideration, shown in Fig. 1, is a large rectangular tank (246 125 720 mm) with a small side branch of square cross-section (6 6 mm) and length 120 mm. Air and water at room temperature are contained in the tank with a flat interface. The height of the interface, h, is measured from the centerline of the branch (h is positive if the interface is located above the branch centerline and negative if the interface is located beneath the branch centerline). The boundary conditions at the top and bottom of the tank, and at the branch outlet For the tests designed to determine the critical height for OLE, a zero pressure was imposed at the branch outlet, a steady pressure P o (whose value depended on the desired Fr) was imposed at the top surface, and water was introduced slowly through the bottom of the tank raising the interface level at a rate of 5 mm/min. Starting with an interface sufficiently lower than the branch, a transient solution was advanced until liquid started to flow into the branch and the value of h at this instant was set as h OLE . Several tests of this type were conducted using various values of P o (each P o corresponds to a different Fr). For the tests designed to determine the critical height for OGE, a zero pressure was imposed at the branch outlet, a steady pressure P o (whose value depended on the desired Fr) was imposed at the top surface, and water was introduced through the bottom of the tank at a rate equal to 99% of the rate exiting through the branch, thus slowly lowering the interface level. Starting with an interface sufficiently higher than the branch, a transient solution was advanced until air started to flow into the branch and the global mass and momentum imbalances were less than 0.01% and the velocity and volume fraction at a selection of points located throughout the branch had levelled off to a constant value.  control volumes using a structured, hexahedral mesh. The mesh was refined towards the interface, the branch inlet, and in the branch itself. ANSYS CFX 14.5 was used in solving the discretized form of the governing mass, momentum, and volume-conservation equations. The k-ε model was used to model the turbulence in both phases and an Eulerian-Eulerian approach was selected treating the flow as an inhomogeneous mixture with each fluid For a fixed value of Po, Fig. 2 shows the dependence of the mass flow rate from the branch on the location of the interface, h. At high interface levels above the branch, the discharge from the branch will be in the form of singlephase liquid. Lowering the interface level, a critical height is reached where the onset of gas entrainment (OGE) occurs at the branch (h = h OGE and m m L,OGE ). A further lowering of the interface results in a two-phase having its own velocity and turbulence fields with discharge through the branch at a rate m TP . As h interaction between the phases via interfacial transfer decreases, m TP decreases while the quality of the twoterms. A parameter that was found to have a significant effect on the onset results is the drag coefficients used in calculating the interphase drag force. This parameter was phase mixture, x, increases. Finally, with further lowering of the interface, another critical condition is reached (onset of liquid entrainment) where liquid stops flowing set at 0.05 throughout this study.
into the branch (h = h OLE and m m G,OLE ). Beyond this Mesh independence tests were conducted in order to assess the accuracy of the numerical results. A sample of these results, corresponding to h = -14.98 mm and P o = 120 kPa are shown in Table 1. These results show that the deviation between the medium and fine meshes is 1.43% point, only single-phase gas flows into the branch.

Onset of liquid entrainment
The onset of liquid entrainment height, h OLE , was in m L,out and 0.11% in m G,out , where m L,out and m G,out are determined at four different values of P o , each the liquid and gas mass flow rates at the outlet of the branch. Based on these results, it was decided to use mesh sizes of 1. 3 [4], and Craya [5]. The following correlation was obtained to fit the present numerical data: A c is the branch cross-sectional area, g is the gravitational acceleration, d is the hydraulic diameter of the branch, ρ G is the air density, and ρ L is the water density. Results were obtained at Fr G,OLE = 3.92, 7.74, 23.56, and 32.83. Figure 3 (a) and (b) shows air volume fraction contours from the transient results for Fr G,OLE = 7.74. The red region in these figures represents 100% air and the blue region represents 100% water. From these results it can be seen that in Fig. 3(a) when h = -10.75 mm, the interface is located far beneath the branch inlet and by raising the interface slightly to h = -10.70 mm, Fig. 3(b) shows that the water beneath the branch inlet climbs the wall and small amounts of water exits through the branch (OLE). It is interesting to note that this phenomenon occurs very quickly (within a 0.05-mm rise of the interface), consistent with previous experimental observations [1][2][3][4].
The discrepancy seen in Fig. 4 between the present correlation and the previous correlations is small. In view of the difference in the branch geometry (square in this study and circular in previous studies), the result seen in Fig. 4 is satisfactory.   At the onset, the air volume fraction at this location changed from 0 to 1. The phenomenon shown in Fig. 5 occurred almost instantaneously while the interface level dropped by less than 0.05 mm. As the interface drops close to the onset height, the interface begins to dip near the wall above the branch inlet and a gas cone forms extending from the interface to the branch inlet without contacting the wall. A very slight further drop in the interface height causes the gas cone to break and the air to flow freely along the wall above the branch as it enters the branch inlet. These observations are again consistent with the observations from previous experiments.

Two-phase flow
When the interface is located above h OLE and below h OGE , both air and water flow through the branch. Under these conditions, results for the total mass flow rate, m TP , and mixture quality, x, at the branch outlet were obtained for various interface heights at two tank pressures: P o = 35 kPa and 120 kPa. well, h OGE and h OLE increase as P o increases. Figure 6 shows the numerical results of h OGE /d vs. Fr L,OGE along with the correlations of Micaelli and Memponteil [2], Yonomoto and Tasaka [3] and Hassan et al. [4]. The numerical data were fit by the following correlation, which is shown in Fig. 6   Hassan et al. [4] found that with their experimental data for discharge flow through a single horizontal branch, the effect of P o was absorbed by plotting their h OGE d 0.636 Fr L,OGE .
( 5 ) data in terms of the dimensionless parameters M and H, where From Fig. 6 it can be seen that the present results fall within the predictions of existing empirical correlations. and Based on their experimental data, they developed empirical correlations for M as a function of H and x as a function of H. The present numerical results for P o = 35 kPa and 120 kPa were plotted in Fig. 8 in terms of M vs. H and in Fig. 9 in terms of x vs. H. The correlations of Hassan et al. [4] are shown in these figures as solid lines. From these figures it can be seen that the present numerical data also collapse when plotted using the dimensionless variables defined by Eqs. (6) and (7). In addition, the present numerical results show very good agreement with Hassan et al.'s correlations.

Conclusion
A numerical study was conducted on the discharge from a stratified two-phase region through a small side branch using a commercial CFD code (ANSYS CFX 14.5). Results were obtained for the flow phenomena of OLE and OGE, the critical heights h OLE and h OGE , the mass flow rate m TP , and quality x of the two-phase discharge. These results agreed well in magnitude and trend with existing experimental data. Thus, it can be concluded that ANSYS CFX 14.5 is capable of predicting the flow phenomena for the geometry under consideration. This study is only a first step and it should be followed by more work involving more elaborate geometries. One drawback is that, even for the present simple geometry, the computation time required for obtaining accurate converged results was found to be excessive (in the order of weeks for one data point using 8 computer cores). In order to avoid divergence, it was found that very small time steps (10 -4 to 10 -3 seconds) and large number of nodes (1.3 to 1.8 million nodes) were required. This can pose a problem in industry when dealing with more sophisticated geometries.