Modelling of Passive Heat Removal Systems: A Review with Reference to the Framatome BWR Reactor KERENA: Part II

: Passive safety systems are an important feature of currently designed and constructed nuclear power plants. They operate independent of external power supply and manual interventions and are solely driven by thermal gradients and gravitational force. This brings up new needs for performance and reliably assessment. This paper provides a review on fundamental approaches to model and analyze the performance of passive heat removal systems exempliﬁed for the passive heat removal chain of the KERENA boiling water reactor concept developed by Framatome. We discuss modeling concepts for one-dimensional system codes such as ATHLET, RELAP and TRACE and furthermore for computational ﬂuid dynamics codes. Part I dealt with numerical and experimental methods for modeling of condensation inside the emergency condenser and on the containment cooling condenser. This second part deals with boiling and two-phase ﬂow instabilities. by an to a high pressure nitrogen bottle. Pressure at the experimental test section is maintained within speciﬁcations by adjusting the pressure in the expansion tank. The experiments were performed at a system pressure of 200 kPa, mass ﬂuxes of 50 to 200 kg/ ( m 2 s ) , and inlet temperatures from ambient to 80 ◦ C. Based on the comparison between experimental results and state-of-the-art analytical approaches are developed. modiﬁed the Chisholm two-phase multiplier correlation put Argonne Laboratory small-channel boiling heat transfer.


Introduction
In part I passive decay heat removal concepts for GEN III(+) nuclear reactors and the modeling of condensation heat transfer in the emergency condenser and in the containment cooling condenser were extensively reviewed and discussed [1]. Especially with the containment cooling condenser, however, not only condensation but also boiling is a relevant thermal hydraulic process that determines the reliability of the system. In general in technical systems with thermal hydraulic circuits, either forced or natural, in which phase change by evaporation takes place, two-phase instabilities may occur and can decisively influence the heat removal and thus the dynamics of these thermal hydraulic processes.
Since these systems have safety-relevant functions such as decay heat removal and shutdown, optimal design and knowledge of the mode of operation is fundamental. Recalculations with integral Subcooled nucleate boiling is initiated, when the saturation temperature of the fluid at the surface is exceeded but the liquid is subcooled. This region is characterized by vapor formation at the heating wall and vapor condensation in the subcooled center of the flow channel. When the fluid becomes saturated, nucleate boiling is fully developed until annular flow starts. In the region of the annular flow, there is a convective heat transfer to the liquid film.
Compared to vertical flow, there are some differences in the flow in horizontal or slightly inclined tubes, because the liquid mainly flows in the lower half of the pipe and the steam mostly flows in the upper one due to gravity. At low flow velocities, where the influence of gravity is large compared to the frictional forces, the two phases are separated (laminar flow, wave flow and slug flow). Due to the asymmetrical distribution of the two phases, the heat transfer coefficient changes over the tube circumference.

ATHLET Subcooled Nucleate Boiling
In ATHLET [17], a modification of the Chen correlation is implemented for calculating the heat transfer during subcooled nucleate boiling [18]. The heat transfer correlation consists of two parts, which consider the microscopic and macroscopic heat transfer mechanism (cf. Equation (1)). For these, the heat transfer coefficient is given as with and where p D is a bundle factor and is defined as a ration of the pitch in meter and the outer diameter of heat conduction volume in meter as well. For Re TP < 4 · 10 5 : For Re TP ≥ 4 · 10 5 : S 0 = 0.133 − 0.825 · 10 −7 Re TP − 4 · 10 5 .
For horizontal tube bundles in cross-flow (γ = 90 • ), no suppression of the microscopic heat transfer h mic could be observed, which results in a factor S = 1. For a tube bundle in a parallel flow (γ = 0 • ) we have S = S0 (cf. Equations (6) and (7)). For all other flows (0 • < γ < 90 • ) S is interpolated with a cosine function between S 0 and 1 according to [17]: Nucleate Boiling The Chen correlation (cf. Equation (9)) is implemented to calculate the heat transfer during nucleate boiling under saturation conditions. Compared with subcooled boiling, the correlation for the macroscopic heat transfer coefficient considers the steam mass quality x (cf. Equation (10)) [18]: The microscopic heat transfer coefficient and the factor F are calculated with the same correlations as for subcooled nucleate boiling (cf. Equations (3) and 4) [18]. Concerning film boiling, ATHLET distinguishes between film boiling at droplet flow and film boiling at inverted annular flow. Starting with film boiling at droplet flow, ATHLET uses a modification of a Dougall-Rohsenow correlation [17,19]: This correlation was developed to predict the heat transfer in a vertical tube or inclined tube by assuming that the flow structure consists of a central liquid core and a thin annular film of vapor on the heated wall [19]. The Groeneveld correlation is also implemented in ATHLET to describe the heat transfer during film boiling at droplet flow [17]. This correlation is applicable in fully developed film boiling within different geometries including tubes, annuli as well as tubes and annuli combined [20]: with The third correlation for calculating the heat transfer during film boiling at droplet flow is a Condie-Bengston IV correlation [17]: This correlation was developed to predict the heat transfer coefficient of tube bundles [21].
Regarding film boiling at inverted annular flow, ATHLET uses a Berenson correlation at first [17]. This correlation neglects the effect of vapor velocity and film thickness on the liquid-vapor boundary behavior during film boiling near the minimum temperature difference [22]: The last correlation concerning boiling, which is used in ATHLET, is the Bromley correlation [17]. This empirical correlation for the heat transfer coefficient is based on forced convection film boiling for an upward flow over a horizontal tube [23]: where 2.1.2. RELAP As already described in [1], RELAP (Reactor Excursion and Leak Analysis Program) was developed to simulate the transient behavior of light water reactor coolant systems during accident Energies 2020, 13, 109 6 of 39 scenarios [24]. Regarding boiling, RELAP distinguishes firstly between different geometries. In this paper the correlations for horizontal single tube and especially for cross flow are considered (cf. Geometry 132 and Geometry 133 in [24]). In RELAP, three different mechanisms concerning boiling inside tubes are distinguished: nucleate boiling, transition boiling and film boiling.

Nucleate Boiling
For saturated boiling, RELAP also uses the Chen correlation for calculating the heat transfer coefficient in a single tube with and without cross flow (cf. Equations (1)-(3)). The difference between the implemented models is the calculation of the factor F and the factor S, which is for RELAP: For subcooled nucleate boiling, a superheated liquid layer next to the hot wall, which is a source of vapor, is assumed [24]. The adaption to subcooled nucleate boiling is provided by a modification of the factor F, which was proposed by Bjornard and Griffith [25] as: Transition Boiling For transition boiling a correlation based on the Chen correlation is used. The implemented model considers wall-to-liquid as well as wall-to-vapor/gas heat transfer and finally the heat flux [24]: q tot,w =q CHF · exp −k · (T w − T sat ) 0.5 · M l + 0.0185 · Re 0.83 Pr 1/3 · (T w − T v ) · (21) · 1 − exp −k · (T w − T sat ) 0.5 · M l , with k = max (k 1 , k 2 ) , Film Boiling Film boiling in RELAP is separated into three mechanisms: conduction, convection and radiation. For film boiling conduction, a correlation by Bromley is implemented in RELAP [24]. The correlation describes the laminar conductive heat flux from a horizontal tube to a fluid at rest: Subcooled film boiling conduction is considered by including a modification by Sudo and Murao [26]: The heat transfer coefficient at convective film boiling is based on a Dittus-Boelter correlation for turbulent forced convection [24]: Regarding film boiling radiation, three different interfaces for the heat transfer are described based on Sun et al. [27]: wall-to-liquid (wl), wall-to-vapor (wv) and vapor-to-liquid (vl). The heat flux for each interface is calculated by the following equations: with The emissivities are given as: where L m is a mean path length, and a v and a l are vapor/gas and liquid absorption coefficients, which are defined as: The vapor/gas absorption coefficient a v and the emissivity w of a Zircaloy-wall are taken directly from references for a fixed temperature [ part was retained to minimize potential unphysical oscillations, which may come up for low-pressure conditions [29]. For pool boiling, various correlations such as Steiner and Taborek [30], Cooper [31] and Gorenflo [32] were investigated under saturated conditions at pressures of 2.7 bar (typical for reflooding and passive cooling conditions) and 70 bar (typical for PWR SBLOCA pressure and BWR operating conditions). The correlations were compared to the McAdams correlation (low pressure) and the Levy correlation (high pressure) as recommended by Holman [33]. The Gorenflo correlation [32] was selected because of its satisfying values for low pressure and high pressure conditions as well as the simplicity of the correlation: where Represent the relative effect of the boiling pressure, the heat flux density and the heating surface properties Ra. The correlation relates to a mean heat flux density of qq 0 = 20,000 W/m 2 and the heat transfer condition of h 0 = 5600 W/m 2 K at this heat flux in case of water. Finally the pressure dependency is included by the exponent N (p * ) and the function N (p * ), where p * is the normalized saturation pressure. The effect on the surface parameter is normalized by the cooper surface parameter R a0 = 0.4 µm.

Transition Boiling
The transition boiling regime provides the transition between the wet wall heat transfer of the nucleate boiling regime to the dry wall heat transfer of the film boiling regime. For TRACE, the minimum film boiling temperature is used as the regime transition criterion between film boiling and transition boiling. The correlation at this point is simply the film boiling heat flux calculated for the wall superheat, where the film boiling coefficient can be either that from dispersed flow film boiling or inverted annular film boiling.

Film Boiling
Also in TRACE, three different film boiling regimes are distinguished: inverted annular film boiling, dispersed flow film boiling and inverted slug film boiling. In the case of inverted annular film boiling, the hot surface is separated from the subcooled liquid core only by a thin film of vapor. The void fraction in this regime is less than 0.6. Three different correlations for the heat transfer coefficient are implemented: wall-to-vapor (wv), wall-to-liquid (wl) and wall-to-interface (wi) [29].

of 39
It is assumed that the dispersed flow film boiling is present at void fractions greater than 90%. This regime consists of two parts: convection and radiation. For laminar forced convection, the constant Nusselt number for fully developed flow with a constant heat flux boundary condition is applied [29]: For turbulent forced convection, as recommended by both Bhatti and Shah and Incropera and De Witt, the Gnielinski correlation is selected for implementation [29]: The correlation by Sun, Gonzalez-Santalo and Tien, which is has already been introduced for RELAP, is also implemented in TRACE [29].
To ensure a smooth transition of the heat transfer coefficient between inverted annular and dispersed flow film boiling, TRACE uses an inverted slug film regime to consider the subcooled liquid flowing past the quench front. Therefore, the heat transfer coefficient comprises of the heat transfer coefficients for inverted annular film boiling (IA) and dispersed flow film boiling (DF) [29]: where:

Boiling Heat Transfer Models in CFD
In addition to one-dimensional codes, more and more CFD codes are used in order to simulate the flow behavior during the boiling, as they provide complete flow fields (i.e., velocity, pressure, temperature, etc.) at all locations of the area of interest. To do this, wall-boiling sub-models must be used to account for the heat input and phase change on a heated wall. In the following, a review of works have been done on boiling with different multiphase flow approaches such as the one-fluid volume-of-fluid model (VoF), level set, and the Euler-Euler two-fluid approach is given.
In recent decades, several CFD models were developed to investigate the bubble dynamics based on simulation of a single-bubble. Lee and Nydahl [34] simulated a single-bubble generation during nucleate boiling by use of moving mesh and the generalized arbitrary Lagrangian-Eulerian (ALE) approach. However, the assumption of constant wall temperature and the lack of suitable modeling of the detachment led to an only partially accurate prediction of the bubble dynamics.
Son et al. [35] modeled bubble dynamics with the level set method for tracking the liquid-vapor interface. The calculated bubble dynamics were in good agreement with experimental results. However, the model does not predict the bubble waiting time and frequency because a constant temperature is assumed for the heat transfer surface.
Tryggvason and Lu [36] used direct numerical simulation (DNS) for modeling nucleate boiling and bubble generation. The micro region was not considered and the heat release rate at the phase boundary was defined as: where ε(n) denotes a delta function of the normal coordinate. Figure 2 shows a developing bubble shape and the temperature profile during nucleation boiling on a vertical wall. The hot liquid layer is generated near the hot wall, which leads to bubble growth. Since the flow moves upward, the generated bubble diverts the thermal boundary layer near the wall so that the cold liquid gets near the hot wall. Sato and Ničeno [37] developed a model based on a color/density function that is able to simulate the dry spot underneath the bubble. Moreover, the bubble growth rate and the temperature distribution over the heat transfer surface were predicted and validated with experimental data. Due to the large computational time, the model is limited to a few nucleation sites and only available for nucleate pool boiling. Moreover, the simulation domain is strictly limited to a millimeter to centimeter size of the global boiling process. Later on, Sato and Niceno [38] proposed a new micro-layer model for nucleate pool boiling using an interface tracking method. They defined the micro-layer thickness as a variable which decreases during the boiling and is able to get vanished at the end of the evaporation process. The model was validated with experimental data of Duan et al. [39] and Yabuki and Nakabeppu [40] and showed a good prediction of the bubble growth rate and the temperature field at the interface. By use of the developed models, Sato and Niceno [41] conducted a series of pool boiling simulations covering the range of the nucleate boiling to the film boiling regarding critical heat flux (CHF). Figure 3 shows the simulation results for two heat flux values. For the case with the low heat flux value q = 50 kW/m 2 , isolated bubbles are formed and superheated liquid is detected. By enhancement of the heat flux, the dry spot area increases and more part of the heat transfer surface is covered with the gas phase. However, direct numerical simulations could provide essential information for deriving a correlation that is valid over a wide range of parameters. For industrial applications, simplification of the sub-process models based on microphysics is still required in numerical simulation.  Fig. 8. The values are calculated in the same way as employed in the experiment: i.e. from the temperature at the horizontal (z-constant) plane 3.7 mm below the heat-transfer surface, where a thermocouple had been installed; Tw is then estimated as follows: where ZTC ¼ À3:7 mm is the position of the thermocouple, and Tz¼Z TC is the space-averaged temperature over the plane z ¼ ZTC within the confines of the lateral extent of the simulation model. The average wall temperature approaches an asymptotic value in all cases, except for q 00 ¼ 1500 kW/m 2 . The computation for q 00 ¼ 1500 kW/m 2 was interrupted when the wall temperature exceeded 145°C, the temperature at which the heater burned out in the experiment. Comparisons of the overall heat transfer coefficient measured in the experiment [20] and that calculated are shown in Fig. 9. The wall superheat is defined as T w À T sat . The measured data is scattered; the reason for this remains unclear. In general, the simulation results agree well with measurement, and lie within the scattering of the measurement data. The heat transfer coefficient peaks at q 00 = 1200 kW/m 2 , the so-called CHF.
In the present simulations, a database trendline of the active nucleation site density as a function of the wall temperature is necessary in prescribing the nucleation-site activation temperatures. The trendline measured by Gaertner [20] is available up to T w = 115.2°C for nucleate boiling at atmospheric pressure (Fig. 1b). The highest value of T w computed in the simulations undertaken here was 145°C, this for q 00 = 1500 kW/m 2 . As a consequence, the extrapolation of the measured tabled results was admittedly used in the simulation.

Bubble shape and temperature distribution
The bubble shapes and the temperature distributions at pseudo-steady-state are presented in Fig. 10 for the cases q 00 = 50, 100 and 300 kW/m 2 . The left column shows the bubble shape in the entire domain and T solid distribution; the middle column is

(b)
fluxes are displayed in Fig. 8. The values are calculated in the same way as employed in the experiment: i.e. from the temperature at the horizontal (z-constant) plane 3.7 mm below the heat-transfer surface, where a thermocouple had been installed; Tw is then estimated as follows: where ZTC ¼ À3:7 mm is the position of the thermocouple, and Tz¼Z TC is the space-averaged temperature over the plane z ¼ ZTC within the confines of the lateral extent of the simulation model. The average wall temperature approaches an asymptotic value in all cases, except for q 00 ¼ 1500 kW/m 2 . The computation for q 00 ¼ 1500 kW/m 2 was interrupted when the wall temperature exceeded 145°C, the temperature at which the heater burned out in the experiment. Comparisons of the overall heat transfer coefficient measured in the experiment [20] and that calculated are shown in Fig. 9. The wall superheat is defined as T w À T sat . The measured data is scattered; the reason for this remains unclear. In general, the simulation results agree well with measurement, and lie within the scattering of the measurement data. The heat transfer coefficient peaks at q 00 = 1200 kW/m 2 , the so-called CHF.
In the present simulations, a database trendline of the active nucleation site density as a function of the wall temperature is necessary in prescribing the nucleation-site activation temperatures. The trendline measured by Gaertner [20] is available up to T w = 115.2°C for nucleate boiling at atmospheric pressure (Fig. 1b). The highest value of T w computed in the simulations undertaken here was 145°C, this for q 00 = 1500 kW/m 2 . As a consequence, the extrapolation of the measured tabled results was admittedly used in the simulation.

Bubble shape and temperature distribution
The bubble shapes and the temperature distributions at pseudo-steady-state are presented in Fig. 10 for the cases q 00 = 50, 100 and 300 kW/m 2 . The left column shows the bubble shape in the entire domain and T solid distribution; the middle column is

(c)
fluxes are displayed in Fig. 8. The values are calculated in the same way as employed in the experiment: i.e. from the temperature at the horizontal (z-constant) plane 3.7 mm below the heat-transfer surface, where a thermocouple had been installed; T w is then estimated as follows: where Z TC ¼ À3:7 mm is the position of the thermocouple, and T z¼Z TC is the space-averaged temperature over the plane z ¼ Z TC within the confines of the lateral extent of the simulation model. The average wall temperature approaches an asymptotic value in all cases, except for q 00 ¼ 1500 kW/m 2 . The computation for q 00 ¼ 1500 kW/m 2 was interrupted when the wall temperature exceeded 145°C, the temperature at which the heater burned out in the experiment.
Comparisons of the overall heat transfer coefficient measured in the experiment [20] and that calculated are shown in Fig. 9. The wall superheat is defined as T w À T sat . The measured data is scattered; the reason for this remains unclear. In general, the simulation results agree well with measurement, and lie within the scattering of the measurement data. The heat transfer coefficient peaks at q 00 = 1200 kW/m 2 , the so-called CHF.
In the present simulations, a database trendline of the active nucleation site density as a function of the wall temperature is necessary in prescribing the nucleation-site activation temperatures. The trendline measured by Gaertner [20] is available up to T w = 115.2°C for nucleate boiling at atmospheric pressure (Fig. 1b). The highest value of T w computed in the simulations undertaken here was 145°C, this for q 00 = 1500 kW/m 2 . As a consequence, the extrapolation of the measured tabled results was admittedly used in the simulation.

Bubble shape and temperature distribution
The bubble shapes and the temperature distributions at pseudo-steady-state are presented in Fig. 10 for the cases q 00 = 50, 100 and 300 kW/m 2 . The left column shows the bubble shape in the entire domain and T solid distribution; the middle column is simulation results agree well with measurement, and lie within the scattering of the measurement data. The heat transfer coefficient peaks at q 00 = 1200 kW/m 2 , the so-called CHF.
In the present simulations, a database trendline of the active nucleation site density as a function of the wall temperature is necessary in prescribing the nucleation-site activation temperatures. The trendline measured by Gaertner [20] is available up to T w = 115.2°C for nucleate boiling at atmospheric pressure (Fig. 1b). The highest value of T w computed in the simulations undertaken here was 145°C, this for q 00 = 1500 kW/m 2 . As a consequence, the extrapolation of the measured tabled results was admittedly used in the simulation.

Bubble shape and temperature distribution
The bubble shapes and the temperature distributions at pseudo-steady-state are presented in Fig. 10 for the cases q 00 = 50, 100 and 300 kW/m 2 . The left column shows the bubble shape in the entire domain and T solid distribution; the middle column is  In the studies presented above, CFD simulates the boiling by resolving the interface between liquid and vapor. This requires a large computational effort, which limits the scope of applications for this CFD-approach. Therefore, the Euler-Euler two-phase flow approach has been considered useful for the large-scale application [42][43][44], because it does not resolve the interface, but applies submodels to describe the conservation equations of mass, momentum, and energy. Two challenges exist for modeling the boiling process with this approach. The first is the wall boiling model, which takes the heat transfer to the liquid and the heat required for the vapor generation into account. The second refers to the phase change that occurs in the bulk fluid, such as condensation and evaporation.
Regarding to the first, most of the Euler-Euler CFD models follow the heat flux partitioning approach initially proposed by Judd and Hwang [45] and further developed by Kurul and Podowski [46,47]. In this approach, a given total heat flux Q tot is split into various terms according to a microscopic model concept (see Figure 4). The total wall heat flux Q tot is determined as a sum of three terms: where Q C , Q Q and Q E denote the heat flux components due to single phase turbulence convection, quenching and evaporation, respectively. The terms Q C , Q Q and Q E are modeled as functions of local flow parameters and local wall temperatures. Several investigations have been carried out to determine the heat flux terms such as Tolubinsky and Kostanchuk [48], Kocamustafaogullari [49], Ünal [50], Klausner et al. [51], Zeng et al. [52,53]. However, until 2013 a critical review of the different correlations applied in the heat flux partitioning model has shown that the parameters are still not suitable for a wide range of application for different fluids or pressure levels. These parameters have to be carefully calibrated for the intended purpose [54]. Colombo and Fairweather [55] as well as Ding et al. [56] used a multiscale bubble growth model and considered the force balance equations for predicting bubble dynamics. However, it has been found that the complex bubble model cannot be easily implemented in CFD codes. Ding et al. [57] succeeded in implementing these bubble models in the CFD approach using a library-interpolation method. In the implementation, libraries of bubble parameters were generated under different conditions, which in turn were implemented in the CFD code and interpolated by the local flow parameters. However, the complexity of this method still limits its wide application.
g Q E is consumed by the vapour generation. After the release of the initial bubble, cooler bulk liquid comes in contact to the surface where the bubble was previously attached. This new liquid adjacent to the wall leads to cooling Q Q . This mechanism, which is obviously not present in single-phase flows, is termed quenching. After a waiting time t w , the thermal layer near the wall is reformed, and a new bubble occurs at the same place. This so-called ebullition cycle extends over a time τ = t g + t w . Consequently, the departure frequency f is calculated as f = 1/τ. On parts of the wall, where no bubbles reside, heat Q C flows directly to the subcooled liquid in the same way as in single-phase flow. The complete cycle is presented in Fig. 1. Judd and Hwang (1976) first proposed the heat flux partitioning concept which was further developed by Podowski (1990, 1991). Detailed informations on heat flux partitioning are also given by Yeoh et al. (2008). A given overall het flux Q tot is divided into various components according to a microscopic model concept (see Fig. 1).
Accordingly, the given external heat flux Q tot , applied to the heated wall, is written as a sum of three parts: In the CCC, film boiling may occur that leads to a vapor film covering the wall surface. In this case, the heat transfer to the vapor phase has to be formulated. Lifante et al. [59] extended the heat flux partitioning model according to Equation (52) with the heat flux which is transferred to the vapor phase (Q V ) due to convection: They applied a switch function f (α L ) to determine the heat transfer to liquid and vapor phase by assuming the critical liquid volume (α L,crit ) equal to 0.2: where α L is the liquid volume fraction,h cv is the vapor heat transfer coefficient, T w and T v represent the wall and the vapor film temperature. The calculated results were validated with two experimental data sets of Bartolomej [60] and Hoyer [61] and gave satisfactory results. Bruder and Sattelmayer [62] investigated the extended wall boiling model with experiments performed at the Technische Universität München in Germany. The work was carried out for a one side heated cubic channel at atmospheric pressure. They found that the void fraction at the critical heat flux (CHF) varies between 0. 35 [64]. Between these bubble size groups, mass transfer is assumed both due to bubble coalescence and fragmentation, as well as condensation and evaporation. Krepper et al. [65] developed the inhomogeneous MUSIG model based on the original MUSIG model with the sign of buoyancy changing due to bubble size.
During boiling inside the tubes of the containment cooling condenser, flow pattern transitions may occur due to the increase in void fraction.
In order to simulate different flow morphologies and the transition between them, the interface has to be resolved. For this purpose Hänsch et al.
[66] developed the GENTOP concept based on the inhomogeneous MUSIG approach by adding a group for continuous gas phase (cf. Figure 5). The combination of the wall boiling model and the GENTOP concept is able to detect the flow pattern transitions [67].
introduced by using a three-field twois represented by a continuous liquid s phase dg and a continuous gas phase ed gas dg is modelled in the framework ltiple Size Group (MUSIG)-approach to e size groups and associated velocity persed gas phases, mass transfers beze groups due to bubble coalescence o account by appropriate models. The odel offers the possibility to add more ith own velocity fields in the future.

SIG-model
eous MUSIG-framework the gaseous into a number of N velocity groups rized by an own velocity field. This ene momentum exchange on the bubble nce on the drag and non-drag forces al., 2007). The overall bubble size disnumber of size fractions M j assigned to . Each size fraction i(i = 1, . . . , M) repre- An occupation number f i = a i / f each size group to the total void fracopulation balance model is applied to to consider bubble coalescence and een different bubble size groups. The ly occur inside the size fractions but phases. Governing equations using this y and momentum equations for a gasuming the same density for all gaseous j;VM þM j;TD ð3Þ rfacial forces acting between gaseous from independent physical effects such ce, wall lubrication, virtual mass and mass transfer between dispersed MUSIG-phase and continuous phase has to be modelled dependent on the flow situation. As a preliminary step the new continuous gas phase here is supposed to be the last velocity group included in the MUSIG-framework as demonstrated in Fig. 2. This continuous velocity group consists of only one size fraction and should represent all gas structures that exceed a maximum value of equivalent spherical bubble diameter d dg,max . Mass transfers between the continuous gas phase to every bubble size group of the polydispersed gas phase are possible by including appropriate models. The completion of the GENTOP-concept to a four-field model adding a polydispersed droplet phase is possible in general. The threshold diameter d dg,max has to be chosen according to the grid size using d dg,max = zDx with 1 0 z 0 10. Considering the gap in the interfacial scales mentioned by Tomiyama et al. (2006) the size d dg,max has to be a compromise between the limitations of the averaged MUSIG-modelling and a useful representation of a resolved gas structure. A resolution of the gas-liquid interface is only possible if the continuous gas structures contain several grid cells. When using the Eulerian approach the interface tends naturally to smear over 2-3 cells. For that reason z % 4 was chosen keeping in mind that this parameter involves uncertainties and has to be checked carefully. As the bubble size d dg,max is several times larger than grid size, at this point it is more convenient to finally resolve the gas structures then using averaged models that have been developed for d 0 Dx. Therefore, the switch between dispersed and continuous gas phase is grid dependent. A proper cell size has to be chosen depending on the simulated bubble size distribution and the size of gas structures for Consequently, the Euler-Euler CFD approach has a high potential to simulate the boiling process in the CCC because it can represent nucleate boiling and film boiling as well as morphologies transition almost realistically.

Experiments
In 1995, Gupta et al. [68] carried out experimental investigations to determine the convective boiling heat transfer coefficient of the local flow in small tube bundles. There was distilled water in the tubes with low cross-flow velocities at atmospheric pressure. The schematic diagram of test facility is represented in Figure 6.
The experimental set-up includes a vessel, heaters, a preheater, a receiver unit, a condensing and cooling water system, and the necessary measuring devices. The horizontal tubes are arranged in a vertical array through a large channel in test vessel. The experiments were performed in three distinct sets of bundle arrangements: (1) a single tube in a channel to represent a basic building block; (2) two horizontal tubes placed one above the other at different pitch distances, and (3) three horizontal tubes one above the other at a constant pitch distance. All conducted experiments are carried out under a system pressure of 1.0019 bar. The heat flux and mass flux ranged from 10 to 40 kW/m 2 and 0 to 10 kg/(m 2 s). The effects of the impressed heat flux, cross-flow velocity and tube geometry on the heat transfer characteristics have been investigated. A Chen-type relation has been used to correlate the data on local forced convective heat transfer coefficients of upper tubes with reasonably accuracy. the bundle) in the vessel and the top of the test vessel was connected to the condenser through a pipe. A pressure gauge and a mercury thermometer fitted at the top of the vessel were used to monitor the system pressure and bulk fluid temperature, respectively. The vessel was fitted with a water level indicator and pockets to insert thermocouples to measure the bulk fluid temperature. A 1 h.p. (0.736 kW) pump was used to circulate liquid from the receiver unit to the test vessel, and the rate of circulation was regulated through a gate valve and a by-pass line ; this rate was measured by means of a rotameter. The fluid entering the test vessel was saturated and its temperature was controlled by a heating element provided in the preheater. Power to this heating element was regulated with a variable transformer. The system components and piping were well insulated to prevent heat losses to the surroundings.
A stainless steel (AISI 304) tube of commercial finish having a 19.05 mm o.d. and 0.89 mm wall thickness formed the test heater, the details of which are shown in Fig. 2. Each tube was 245 mm long and, with an effective heating length of 190 mm, provided resistive heating by means of a high alternating current fed to it through a step-down transformer. All tubes of the bundle were connected in series and the amount of current passing through the bundle was controlled through a variable-voltage transformer. Power supplied to test heaters was measured with the help of a digital voltmeter and digital ammeter with a lowest Exiting the flowmeters, the fluid flows through a system preheater in which the fluid temperature is raised to a desired subcooled level for a given test. The preheater consists of a Type 304 stainless steel tube with 4.57-mm inside diameter, 6.10-mm outside diameter, and 500-mm length. Direct current at low voltage from a regulated power supply is passed through the preheater tube wall to generate resistance heat. The power supply output regulation range is 0.1% of the voltage or current and has a power capability of up to 10 kW.
After passing through the preheater, the fluid enters the experimental test section, shown schematically in Fig. 2. The test section is resistance heated with DC by a separate regulated power supply (maximum output: 18 kW) from the preheater. Voltage drop across the test section is measured directly, and the current is determined from a measurement of the voltage drop across a shunt resistor. Power (heat) input to the test section is calculated as the product of the voltage drop and the current. Immediately beyond the experimental test section, a Pyrex window provides a view of the flow pattern. Electrical isolation of the preheater/test-section arrangement, for the purpose of eliminating ground loops, is provided by short sections of high-pressure hose, designated ISO (for isolation) in Fig. 1.
The test section is well insulated thermally from the atmosphere to minimize heat loss to the environment. However, test-section heat loss calibration tests were performed, and the slight heat Wojtan et al. [70,71] built up a test facility to acquire the experimental data of flow boiling in horizontal tubes. Based on the dynamic void fraction measurements described in [72], he developed a new flow pattern map to divide the stratified-wavy region into three subzones: slug, transition regime and stratified-wavy. And the annular-to-dryout and dryout-to-mist flow transition curves were added and integrated as well. Additionally, Wojtan developed a new heat transfer model for stratified-wavy, dryout and mist flow regimes. Figure 8 shows a simplified layout of intube refrigerant test loop with a close up view of the set-up used for the dynamic void fraction measurements.
The refrigerant first goes through a series of horizontal electrical preheaters and then passes an insulated tube without any sharp elbows. Then, the refrigerant enters the tubular test section and is heated by counter current flow of hot water in the shell side of the double pipe system. The temperature of the hot water is measured at five positions (using 26 thermocouples) that allows the enthalpy profile along the test section to be determined and from that the local heat flux q. The temperature of the refrigerant T sat is determined by the measurement of the saturation pressure at the inlet and outlet of the test section (assuming a linear variation). The wall temperature T wall is measured using four thermocouples installed on the external wall of the tested tube (taking the mean of the measurements at the top, two sides and bottom of the tube wall). Then, the refrigerant local heat transfer coefficient was calculated as follows: where D ext and D are the external and internal diameters of the copper tube and k is the thermal conductivity of copper.
More technical details concerning the test section and the heat transfer measurements are presented in Part II of this paper.

Development of the new version of the flow pattern map
The flow pattern oriented heat transfer model of Kattan et al. [4] was developed for vapor qualities higher than 0.15. According to the flow pattern map presented in Fig. 6, all points less than the vapor quality of 0.15 and below the G wavy curve are identified as stratifiedwavy and for the lowest mass velocities (below G strat ) as stratified flow. Kattan [3] suggested that the frequency of waves in this region is very influential on the local heat transfer coefficient. He proposed in the future to divide this zone into two subzones according to an all wet or partially wet criterion. Zü rcher [16] found also that the flow pattern oriented model of Kattan et al. did not work as well for stratified-wavy flows as for annular ones.
In this light the dynamic void fraction results from [2] will be used to discern the transitions between the respective flow regimes in stratified types of flow. The optical void fraction measurements technique is applicable only for stratified, stratified-wavy and slug flows. In the present tests, the slug and wave slopes were nearly never higher than the inclination angle of the camera, which would block the cameraÕs view of the interface. It allowed using the optical transformation developed for the static conditions taking into account only light  The schematic diagram of the experimental facility is shown in Figure 9. The test facility includes a semi-hermetic compressor, a water-cooled condenser, an expansion device, a preheater, a postheater and a evaporator. The experiments were carried out at inlet temperatures from 6 to 9 • C, heat fluxes between 3 and 6 kW/m 2 and refrigerant mass fluxes from 100 to 300 kg/(m 2 s). During the boiling of the refrigerant steam mass qualities between 0.1 and 0.9 were reached. Based on the experimental results an empirical correlation for the prediction of the heat transfer coefficient of R407C during boiling inside inclined tubes was developed. rate has been controlled through hand shut off valves fitted on the by-pass line incorporated after the condenser. A quality filter-dryer was accommodated after the flow meter to entrap lubricating oil, foreign particles and moisture in refrigerant. A suitable pre-evaporator was designed and installed to control the vapor quality at the inlet of test evaporator. By regulating power supply with variable voltage AC, heat source, heat input to pre-evaporator has been controlled. An oil separator was incorporated downstream of the compressor for oil free refrigerant to enter in the test evaporator section. A suitable accumulator was also installed upstream of the compressor to ensure that the oil free refrigerant vapor should enter the compressor.
The test section was made of smooth copper tube with inner diameter of 7 mm and outside diameter of 9.52 mm. The outside tube wall temperatures were measured by T-type copper-constantan thermocouples at five axial stations, including inlet and exit to the test evaporator tube. At each station, the temperatures were measured at the top, two sides of the middle and bottom position. The average of these temperatures indicates the local wall temperatures. The local saturation pressures at the inlet and outlet of test evaporator were measured by piezoelectric pressure transducers having a range of 0e20 bar. The average of these two pressures indicates the saturation pressure at test section. The test evaporator tube was heated by flexible Nichrome (Nickel-eChromium 80e20 percent by weight) heater wire of 3.2 kW capacity at 8 m length wrapped over mica thin tape which was wound around the tube over the full test length of 1.2 m. Heat

Data redu
The 'quasi-local' h using following Equ This is simply resistance heating and conduction h refrigerant flowing outside tube wall te temperature at the the tube at each se T Sat is the satu pressure readings, point temperature the mean of the lo the heat loss to t defined as the rat power, the resista outside of the tub calculated by using x (typically arou using the method Average vapor qua where, x in is the test section, which

Two-Phase Flow Instabilities
For the passive removal of the decay heat from the containment into the storage pool above by the CCC, the principle of natural circulation is applied. Depending on the supplied heat, the fluid passes from a single-phase flow into a two-phase flow by the boiling process. The phase change from single-phase flow to two-phase flow initiates two-phase flow instabilities, whose phenomenology is known in both natural and forced circulation systems. The instability mechanisms are associated with the propagation of perturbation in a two-phase flow due to feedback effects, which involve delays [74]. Generally, these mechanisms are divided in two mechanisms; first, microscopic mechanisms, which address instabilities occurring at liquid-gas interfaces and are neglected in this article, and second, macroscopic mechanisms, which refer to instabilities in the whole system. Bouré published the most cited classification of the macroscopic instability, the static and dynamic mechanisms as first distinction [8] (c.f. Tables 1 and 2). The Sections 3.3 and 3.4 address a brief overview of experimental and numerical results concerning two-phase flow instabilities in natural circulation systems. The static instabilities mainly differs from the dynamic instabilities, that the thermal dynamic equilibrium of the two-phase flow is restored after perturbation in case of static instabilities. Concerning static mechanisms, the unstable threshold can be predicted by solving the steady state conservation equations. Whereas the various dynamic feedback effects which occur regarding dynamic mechanisms influence the threshold of unstable behavior significantly.

Ledinegg Instability
The Ledinegg instability is well known and intensively investigated. Ledinegg introduced this type at first in Ledinegg [75].
The flow excursions occur when the slope of the channel characteristic for boiling systems is more negative than the external characteristic at small perturbation around an equilibrium point [8]: Figure 10 shows the trend of the internal pressure drop in the cases of all vapor (x = 1), all liquid (x = 0) and mixture (0 < x < 1) in a cross-flow channel. Furthermore, two cases of external pressure drop trends induced by several pumps are shown. In case one, the characteristic curve of the pump has three intersections within the mixture characteristic curve. The operating point B is unstable according to the stability criterion (cf. Equation (59)) and a small perturbation causes a change to the operating point A or C. If the new stable operating point A is reached, the system will be at the risk that the necessary cooling of the heated wall is not guaranteed and the critical heat flux is reached. This instability can be avoided by installing an inlet throttle valve, which leads to the characteristic curve of case two [8]. The stability criterion according to Equation (59) is now fulfilled and operating point B is stable. The external pressure drop of a natural circulation system is the buoyancy due to the gradient of the density. Whereas the internal pressure drop includes all losses in the inlet, boiling channel, exit and downcomer but the pressure drop due to gravity.
Detailed studies of the nonrecurring excursive instabilities with experimental background are published by Maulbetsch and Griffith [78], Bouré et al. [8] and Padki et al. [77]. Padki et al. discovered that the static Ledinegg instability is caused by a saddle-node bifurcation [77].
Nayak et al. [79] studied the stability behavior of the AHWR concept by prediction of the instability threshold in 1998. They concluded, inter alia, that the Ledinegg instability will disappear, if the system pressure is higher than 7 bar or the subcooling less than 10 K during start-up.
Rao et al. [80] carried out a linear stability analysis in a frequency domain in 1995. They figured out that the Ledinegg instability does not occur at the simultaneous presence of a neutronic and hydrodynamic feedback.
Ruspini et al. [81] investigated this phenomenon under low flow conditions by applying model order reduction (c.f. Figure 11). The solution of steady state results in the typical N-shape curve. If the slope is smaller than the external characteristic, the point will be unstable. Two simulations with different perturbation lead to a lower (in the case of an outlet pressure of 996.51 kPa) or higher (in the case of an outlet pressure of 996.49 kPa) mass flow rate.

The Boiling Crisis
The boiling crisis is caused by a non-wetted wall of a boiling channel and is characterized by an insignificant degradation of the heat transfer. At systems where the heat flux is fixed (e.g., nuclear reactors or electrical heating), the non-wetted wall leads to a steep increase of the wall temperature. This can even lead to the failure of heated surface (e.g., burn-out). If the heat is supplied by a given temperature gradient (e.g., steam heating), the heat flux will be strongly decreasing [82]. A review of the boiling crisis in nuclear reactors is given by Theofanous in 1980 Theofanous [83]. There are two existing types of the boiling crisis. The first type occurs by exceeding a certain CHF at a low vapor quality. The steam bubbles coalesce and form a steam film on the surface of a heated wall [84,85]. In case of a higher flow velocity and pressure, the small steam bubbles break away from the surface and form a high viscose bubble layer near the wall. This layer leads to the wetting of the heated wall. The second type of the boiling crisis occurs at the end of an annular flow (c.f. Figure 1). The liquid layer on the wall surface becomes very thin and disappears finally. This causes the degradation of the heat transfer because the heat is transferred to the steam.
Kim et al. [86] investigated the influence of the flow oscillation on the CHF under low power and low flow conditions in 1999. Nikolayev et al. [87] presented experimental investigations of the layer forming and concluded that their results are not compliant to the theoretical aspects of the Zuber model and the macrolayer evaporation theory.

Flow Pattern Transition Instability
The flow pattern transition instability is classified as fundamental relaxation instability. As shown in Figure 1, there are various configurations of gas and liquid flow in a channel, in which these phases arrange themselves. These characteristic configurations are known as flow regimes. Instability occurs during the transition of bubbly to slug or annular flow. A small perturbation or reduction in the volumetric flow increases the quality, which changes the flow regime to the annular flow. The annular flow is characterized by a significantly lower pressure drop than the slug flow. Therefore, the flow rate increases and hence the quality decreases, these two effects lead to the happening of transition again [8].
The great challenge concerning the investigation of this type of instability is the prediction, when the slug flow is initiated. Griffith and Wallis [88] proposed an empirical flow regime map in 1961. Haberstroh and Griffith [89] figured out that the transition depends on the liquid flow rate and created an annular-slug transition criterion. Krussenberg et al. [90] investigated the two-phase flow concerning the transition to slug flow with a reasonable agreement to Taitel flow maps experimentally in 2000 (c.f. [91]). Krussenberg concluded that the slug flow occurs, when the bubble size distribution reaches higher values than the tube diameter. Further, analytical studies are presented by Nayak et al. [79] and Jeng and Pan [92].

Geysering
In general, the compound relaxation instabilities are non-periodical trends of flow dynamics and are irregular. The prediction of these instabilities bases on a steady state solution. General characteristics (e.g., amplitudes, frequencies) cannot be estimated analytically. The geysering instability is a classic example of a compound relaxation instability (c.f. [8]). By reaching the boiling boundary within the boiling channel, the increasing void fraction leads to a decreasing hydrostatic pressure along the channel. This results in an additional boiling above the detached steam bubble with an associated increasing mass flow. The steam bubble is condensing at the inlet of vessel filled with subcooled liquid and causes an eruption in the channel outlet.
Thereupon, the overlying subcooled liquid returns to the channel, which leads to the stagnation of the natural circulation flow [93]. This phenomenon was reported for the first time by Griffith [94]. He observed periods between 10 s and 100 s and finally concluded the occurrence of it at a low pressure and a low power.
Ozawa et al. [95] investigated this phenomenon experimentally on a forced circulation loop for the refrigerant R-113 in 1979. The three main processes of geysering reported by Griffith [94] (single-phase flow, ejection of vapor-liquid mixture and backflow) have also been observed. The period and the amplitude increase with an increasing heat flux and a decreasing inlet velocity. An increasing riser length also results in an increasing geysering-amplitude [95]. Aritomi et al. [96] also discovered in 1990, that the geysering period correlates with the boiling delay time, τ bd , which is the required time for the boiling of subcooled water, which flows through a boiling channel. The boiling delay time is defined as: where ρ l is the liquid density, i sat,l the liquid saturation enthalpy, i l the inlet liquid enthalpy and q the volumetric heat generation rate.
In 2007, Marcel [97] observed only geysering instabilities in the test-facility CIRCUS with a parallel chimney configuration. At low pressure, the phenomenon was observed to be the most prevalent during the investigations on the open natural circulation test facility GENEVA by Cloppenborg et al. [98].
In Figure 12, two-phase-flow oscillations are shown, which are induced by geysering. If the large slug bubbles reach the riser outlet, the pressure drop in the riser will be reduced by the frictional pressure drop and the gravity of about 670 mbar. The gravity pressure drop results from the riser length of about 6 m. The time lag between the mass flow and the riser inlet temperature is nearly 18 s. During the reversal flow, the riser inlet temperature decreases to the subcooled temperature 90 • C and increases during the first phase according to the geysering phenomenon (c.f. Figure 13).

Dynamic Instabilities
Dynamic instabilities were classified according to their causes for propagation of disturbance: density wave oscillations (DWO), pressure drop oscillations (PDO), thermal oscillations and acoustic oscillations. Table 2 represents the classification according to Bouré et al. [8]. Generally, there are two types of transportation of disturbances: pressure and void waves.

Acoustic Oscillation
Acoustic oscillations, also described as pressure wave oscillations, have been observed in the subcooled boiling, the bulk boiling and the film boiling. These oscillations appear at rather highly subcooled operating systems. A high frequency of about 10 Hz to 100 Hz was reported by several investigators [8]. The DWO produces lower frequency waves of about 1 Hz. Bishop et al. [100] detected audible frequency oscillations of 1000 Hz to 10,000 Hz during their experimental investigations in near-critical temperature and supercritical pressure operational conditions in 1964. The amplitudes are relatively small at such high frequencies. The thermal feedback of the vapor film to the passing pressure wave is assumed to be the mechanism for the oscillations during film boiling [12]. Detailed studies are published in [101][102][103].

Density Wave Oscillation
The density wave oscillation is the most common type of flow instabilities. The mechanism is clearly understood and consists of two sub-mechanisms concerning the generation and propagation of the disturbance: the delay of the propagation and the feedback effects on the inlet conditions [104]. Fukuda and Kobori [105] classified these oscillations into two main types. The Type-I (DWO I ) occurs at a very low exit quality, whereas a high quality at the exit is the characteristic feature of Type-II (DWO II ). The pressure of the system is the key indicator for the dominant DWO type. Systems with a pressure below 20 bar (e.g., reactor start-up) are defined as low-pressure system and high pressure systems are specified with a higher pressure than 20 bar [106].
The saturation conditions in a low-pressure system strongly depend on the pressure level, which is shown in Figure 14b,c.
DWO I types are low frequency oscillations in a low pressure systems (c.f. Figure 14b) with a low exit quality. The gravity pressure drop of the adiabatic chimney dominates the system and is very sensitive to flow rate fluctuations. Therefore, the length of the riser plays an important role. With a small perturbation concerning the quality, the void fraction and finally the driving head undergoes a large change. An increasing heat supply suppresses the fluctuation of the driving head for small changes in the quality because of a decreasing slop of the void fraction versus quality [12]. The high frequency DWO II is important for a high pressure (c.f. Figure 14a) and a power, at which forced or natural circulation reactors operate. The oscillations are caused by the interaction of single and two-phase pressure drops, the mass flow and the void fraction in the two-phase region. A different propagation time delay in the single and two-phase region leads to pressure drop variations, which are out-of-phase.
Furthermore, the high-order DWO III is analytically explained by Yadigaroglu and Bergles [107] and Yadigaroglu and Bergles [108]. During the experiments under atmospheric conditions, the DWO III appears at a high subcooling and low power. In 1985, Achard et al. [109] applied the stability island to describe this phenomenon.
Fukuda and Kobori [105] subdivided these types into gravity or frictional effects of a heated section or riser (DWO I-R , DWO I-H , DWO II-R , DWO II-H ) in 1979. Based on their experimental results in [110], these different types of instabilities are obtained by a linear stability analysis in a frequency domain and were observed at different configurations.
An important appearance, where density wave oscillations are induced, is flashing, which is also known as adiabatic boiling. The common flashing induced oscillations appear during a low-pressure operation, for example at the start-up of BWRs (c.f. Figure 14c). The steam is generated in the adiabatic riser because of the decreasing gravitational pressure drop and the saturation temperature along the riser. The steam generation causes an increase in the flow rate and a decrease in the hydrostatic head. The higher flow rate reduces the quality in the heated and riser section resulting in the upwards motion of the boiling boundary. The temperature in the channel becomes lower and the flow rate decreases, which leads to a higher dwell time in the heated section. The process starts again. Flashing is affected by the hydrostatic head and hence sometimes referred as DWO I . The phase difference of the temperature between the heated and the adiabatic section is similar and the dwell time of the fluid and the oscillation period is nearly the same [111]. This phenomenon was experimentally investigated by many researchers such as Schuster [112][113][114], Furuya et al. [5,115], Manera [116,117] and Cloppenborg [98,99].

Thermally Induced Oscillations
The thermally induced instabilities mostly appear in combination with DWO and are the result of the occurring critical heat flux. The high frequency DWO effects a disturbance in film boiling [8]. Due to the changes concerning the flow regimes (c.f. Figure 1) and finally in the heat transfer mechanism, the occurring critical heat flux is moving downward or upward the boiling channel. These

Parallel Channel Instabilities
The parallel channel instabilities can occur in single and two-phase systems, in which the channels are connected to a header. The pressure difference is same for all the channels. In and out of phase oscillations are observed, whereby out of phase oscillations play a dominant role. Kakac et al. [118] found a phase shift of 180 • during experiments with two channels in 1974. In 1964, Berenson discovered a phase shift of 72 • for five channels [119]. Summarizing the experimental results, the general statement of a phase shift with 2π/n is formulated, where n is the number of channels [120].

Pressure Drop Oscillations
In general, pressure drop oscillations are caused by interactions between the channel and a compressible volume (e.g., surge tank or pressurizer) at the inlet of the heated channel and occur commonly at forced circulation system, because of the destabilization by the pump. Contrary to DWO, PDO occurs at higher flow rates and smaller frequencies. PDO is able to stand the same amount of the risks of occurring CHF as Ledinegg.

Natural Circulation Instabilities
According to Aritomi et al. [121], natural circulation oscillations occur during an increasing heat supply, when geysering is suppressed. This in-phase instability is induced by fluctuations of the hydrostatic head in an adiabatic long channel and disappears with a further increasing heat supply and vaporization rate.

Experiments
The experimental investigations provide the basis for the evaluation of the analytical investigations with respect to the modeling. The following section gives a brief overview of the experimental investigations of natural circulation loops concerning two-phase flow instabilities.
Wallis and Heasley [6], the most cited publication, investigated the oscillatory behavior in a natural circulation loop experimentally (cf. Figure 15) and numerically. They discovered DWO I and DWO II during their natural circulation loop operated with pentane. Dijkman et al. [122] investigated the stability of a water operated natural and forced circulation system in 1967. The transfer functions have been measured for the investigation of the stability at steady state conditions. The heated section is electrically heated. The system pressure was varied between 2, 14.7 and 30.4 bar and the power supply was adjusted up to the burnout. The estimated low frequencies of about 0.9 Hz suggest the occurrence of DWO II .
A vast parameter analysis was published by Mathisen [123], who varied the system pressure, inlet and outlet restriction of a single and parallel boiling channel. Mathisen concluded, that the increase of the system pressure, the initial inlet restriction in the separate channels and the ratio of the effective head to the riser length have definitely stabilizing effect regarding the onset of Ledinegg phenomenon. The best performance was observed at minimum subcooling and minimum main flow restriction with parallel boiling channels.
The classification of DWO according to Fukuda and Kobori [105] is based on their experimental investigations. The results were published in [110]. They observed DWO I , where the gravitational factor is the governing factor, and DWO II , where the frictional pressure drop is dominant. It was also shown that the PDO is included in the DWO I .
Schuster investigated the transient behavior of a two-phase natural circulation, which occurs in the test facility DANTON. This test facility represented the axial main dimension of the soviet nuclear district heating station AST-500 in scale 1:1. The start-up was performed in 1985 with respect to the investigation of the thermal hydraulic behavior [113].
In 2000, Schuster et al. [124] stated that the passing through an unstable two-phase region to reach a stable two-phase flow is unavoidable without external pressurizing. During the start-up three types of instabilities are observed: flashing, geysering and DWO II . The main results generally agreed with [96,125]. The study of the start-up procedure was not relevant for the future because the AST-500 is not pursued anymore. But the detailed study of the flashing phenomenon published in [112][113][114] was not yet been done before. The observed flashing at a low pressure near ambient conditions and a low heat flux density induced oscillations with frequencies lower than 0.5 Hz. Figure 16 shows results of the simplified model, which bases on the steady-state equations of energy, momentum and mass. This dependence agrees with the experimental results qualitatively. The riser inlet subcooling has a greater influence on the flashing onset than the riser inlet temperature. Figure 16. Onset of flashing depending on the riser inlet temperature and subcooling [114].
It was shown that this phenomenon has been taken into account for a sufficient start-up procedure and an effective reactor configuration in natural circulation BWR. A detailed discussion of a natural circulation BWR is introduced, inter alia, by Aritomi et al. [96,121]. The experimental investigation of a natural circulation boiling process is carried out in a single vertical boiling channel [96] and parallel boiling channels [121] under the atmospheric pressure. Based on the low-pressure natural circulation  [121] at a heat flux density up to 800 kW/m 2 in a 1 m long heated section, Aritomi et al. estimated three types of two-phase flow instabilities. With an increasing heat flux density, the natural circulation rate is increasing and the geysering is initiated. When the inlet velocity enhances due to an increasing heat flux, the geysering is superseded by a so-called natural circulation oscillation induced by the hydrostatic head fluctuation (cf. Figure 17). The density wave oscillations are following them. At a high subcooling (N sub = 156) and a low power (N pch = 16), a chaotic behavior of the low pressure natural circulation loop was detected by Wu et al. [126]. This results are nearly compliant with [124].
A further natural circulation nuclear reactor is developed by the Institute of Nuclear Energy and Technology (INET) and is in operation since 1989. This 5 MW nuclear heating reactor (NHR-5) operates at system pressure of 15 bar. Jiang et al. [127] investigated the start-up at HRTL-5, which represented the geometric and the system design of NHR-5, experimentally in 1995. They observed geysering, flashing (at a pressure lower than 3 bar) and DWO I during the start-up and proposed a stable start-up procedure. The intense research on the physics of the static and dynamic behavior of the natural circulation operating Dodewaard reactor is published by [106,[128][129][130]. Table 3 shows a brief overview of the experimental investigations in the test facilities with respect to the two-phase flow instabilities in BWR's.
Large integral test facilities such as INKA, PANDA, PUMA and KATHY were established to observe various feedback effects and boundary conditions, which are realistic in BWRs.
The results of the experimental investigations contain the following main statements to stabilized the two-phase flow in general, which are also summarized by [78]: • Increasing the single-phase flow pressure drop in the heated section • Increasing the system pressure • Decreasing the inlet subcooling

Numerical Investigations
The boiling of a fluid flow within a channel is characterized by a two-phase flow and the resulting disequilibrium. The various spatial and time-dependent distributions of the state variables in the vapor and liquid effects stochastic interactions at the phase boundary interfaces. The analytic formulations of the complete complex system require a spatial averaging as a simplification due to their large computational effort. With averaging the state variables for an one-dimensional flow in general, three conservation equations (mass, momentum and energy) for each phase exist. The complexity increases at the transition from a mixture to multi component model. Mixture models are developed for describing a two-phase flow (liquid and vapor) within a defined validation range and assumptions as a pseudo-single phase fluid.
The first most famous cited lucid treatment of the natural circulation instabilities was presented by Wallis and Heasley [6]. They studied the effect of the inlet throttling of the two-phase flow instabilities and concluded that the increasing single-phase pressure drop stabilized the system. A linear stability analysis using a homogenous equilibrium model based on a lumped parameter system was carried out by Zuber [139]. He has given the first introduction of the drift velocities and finally the first approach of the drift-flux mixture model. For example, he discovered that the gravity terms in the single-phase heated section are damping, whereas the gravity terms in the two-phase section are driving terms [125].
Beside the classic method of linear frequency-domain stability analysis, which has been used to study density wave oscillations in [7], the nonlinear stability analysis has attracted considerable interest since 1990s. The nonlinear analysis was used at first by Achard et al. [109], Uddin and Dorning [140] and Clausse and Lahey [141]. The nonlinear effects of two-phase flow dynamics have also been investigated by Pinheiro Rosa and Podowski [142]. The observed results strongly depend on the numerical approaches and the applied spatial discretization, especially for the operating conditions in the linear unstable regions. Pinheiro et al. compared the effects of perturbation on various two-phase flow approaches: the homogenous equilibrium model and two different models of subcooled boiling, which seems to be a significant stable criterion. Furthermore, the results show, that the HEM agrees well with the experimental data at a low subcooling number and phase change number, whereas the drift flux model approximates more exactly at a higher subcooling number and phase change number [142]. Dokhane [143] concluded that the influence of subcooled boiling on the dynamic behavior is sensitive at a low subcooling. Table 4 represents a brief overview of numerical investigations of two-phase flow instabilities in boiling channels or of a whole BWR system. The nonlinear effect was investigated in detail by application of various techniques of model order reduction, especially with the respect to the bifurcation theory. This theory is explained in detail by Hassard et al. [144]. The advanced BWR-ROM of Lange et al. [145] and the application of the nonlinear stability analysis shows a good agreement with the oscillatory behavior of several operating points of the Swiss nuclear power plant in Leibstadt. But the application of the weighted residual method of DFM based on a two-phase model requires a high computational effort. In contrast, the proper orthogonal decomposition (cf. Prill [146]) offers the opportunity to obtain the ROM based on a set of data time series and provides the best approximation.   [154] Linear and nonlinear HEM Reduced order model for the natural circulation and the forced circulation BWR • Low pressure and high pressure systems are modeled • Evaluation with the experimental data of SIRIUS-facility [111] and Dodewaard reactor [106] with nuclear coupling Lange [155] Linear and nonlinear DFM BWR reduced order model with nuclear coupling • Approach of Dokhane [143] extended with recirculation loop model and model feedback reactivity coefficients • Model order reduction with the weighted residual method Ruspini [104] Linear HEM Numerical and experimental investigation of Ledinegg, PDO and DWO • Least squares spectral method is used Prill [146] Linear In 2017, Pandey and Singh [159] published the results of the nonlinear stability analysis of a high pressure natural circulation loop based on previous verified and developed ROM by Paul and Singh [157]. The formulated two-phase region is based on the drift-flux mixture model. The approaches are comparable with those of [150], whereupon they use a linear approximation for the single phase density. Super-and subcritical bifurcation were detected at certain design parameters, especially when DWO I and DWO II stability boundaries diverged from each other [157].
In Figure 18, the DWO are represented by a Hopf bifurcation and the Ledinegg instabilities are characterized by a saddle-node bifurcation [159]. At low-pressure conditions, the pressure dependency of the saturation condition is pronounced (cf. Figure 14b). Zhou introduced a linear approximation of the saturation enthalpy as function of the pressure in the single-phase region [154]. But his two-phase region considered a homogeneous mixture model. This aspect and the high sensitivity of the saturation condition are the keynotes for a realistic formulation of a ROM and finally the application of nonlinear stability analysis.     Fig. 4).

Conclusions
The KERENA TM concept pursues the strategy to simplify the systems, to achieve a higher reliability and a higher safety standard by the application of passive safety system for shutdown and decay heat removal. The simulation of these systems by integral codes, however, shows potential for further development, since the current empirical correlations for the heat transfer at the occurring low mass flows are not reliable. The purpose of this paper (part I and part II) is to summarize and discuss these existing empirical correlations to the relevant thermal hydraulic processes in the emergency condenser (EC) and the containment cooling condenser (CCC). These relevant thermal hydraulic processes can be categorized as: steam condensation inside incline tubes (EC), steam condensation on inclined tubes (CCC), boiling on inclined tubes (EC), boiling inside inclined tubes (CCC) and the natural circulation as a driving force. The first part covers the first three topics [1]. This article is the continuation of the first part with a focus on the boiling process inside inclined tubes and the two-phase instabilities in natural circulation systems. The correlations of the heat transfer implemented in the one-dimensional integral codes ATHLET, RELAP and TRACE are introduced and discussed. Several experimental investigations and the used test facilities are described and summarized, which have formed the base for the formulation of the heat transfer correlation in this special case of low mass flux. In addition, an overview is given to already performed simulations using CFD. Furthermore the scientific work on analytical and experimental investigations of the two-phase instabilities is summarized.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.