Thermal hydraulic considerations of nuclear reactor systems: Past, present and future challenges

Thermal hydraulic analysis of nuclear reactor core and its associated systems can be performed using analysis system, subchannel or computational fluid dynamics (CFD) codes to estimate the different thermal hydraulic safety margins. The safety margins and operating power limits under different conditions of the primary as well as secondary cooling system such as the system pressure, coolant inlet temperature, coolant flow rate, and thermal power and its distributions are considered as key parameters for thermal hydraulic analysis. Considering the complexity of rod bundle geometry, boiling heat transfer and different turbulent scales bring about the many challenges in performing the thermal hydraulic analysis to ensure the safe design and operation of nuclear reactor systems under normal and abnormal conditions. A comprehensive review is presented of past, present and future challenges in state-of-the-art thermal hydraulic analysis c overing various aspects of experimental, analytical and computational approaches.


Introduction
As consideration of the effect of greenhouse gas emission to the environment becomes ever more important, the contribution of global energy requirement by nuclear power is gaining significant traction. Following the nuclear emergencies caused by Three Mile Island (1979), Chernobyl (1986 and Fukushima (2011), a concerted effort to improve the passive core cooling has resulted with advanced nuclear reactor systems with enhanced safety features. Lessons learnt from the accidents of these nuclear systems should be considered in the new design of reactors, especially the reactor core and its associated systems. From thermal hydraulic considerations, the safety requirements of the nuclear reactor systems to be ensured must be able to perform under all different conditions during normal operation, operational transients, anticipated operational occurrences, design basis accidents and under extreme emergency situations by incorporating the engineered safety systems by passive means. This can usually be achieved via thermal hydraulic analysis (Sha, 1980) where such analysis is performed using analysis system, subchannel or computational fluid dynamics (CFD) codes to estimate the various thermal hydraulic safety parameters like critical heat flux (CHF), critical power, fuel centre line temperature, fuel surface temperature, subchannel maximum temperature and bulk coolant outlet temperature (Chelemer et al., 1977). In this paper, the past, present and future challenges of thermal hydraulic analysis for nuclear reactor systems are reviewed based on the design and operation of existing Generation II, III and III+ and advanced Generation IV reactor technologies.

Historical considerations of the predictive capability for nuclear reactor systems
Commercially operated Generation II, III and III+ reactors for power generation are nuclear reactors which they either adopt the Light Water Reactor (LWR), Heavy Water Reactor (HWR) or Pressurised Water Reactor (PWR) technology. From the aspect of performance as well as safety, the understanding of the boiling heat transfer in the pre-critical heat flux region is an important consideration especially with reference to its upper limit. The term Boiling Crisis provides adequate explanation of the phenomena in a symptomatic Vol. 1, No. 1, 2019, 3-27 Experimental andComputational Multiphase Flow https://doi.org/10.1007/s42757-019-0002-5 G. H. Yeoh 4 sense and CHF is commonly associated with the dropping-off of the heat transfer effectiveness in sufficiently cooling and moderating these types of nuclear reactors (Theofanous, 1980).

Flow boiling crisis
To further examine the state-of-the-art in boiling heat transfer particularly in the pre-CHF region, the Flow Boiling Crisis (FBC) can be illustrated in the schematic diagram in Fig. 1 as demonstrated by the simultaneous flow of water and steam that can result in many different forms of the liquid flow and heat transfer as it is transported along the vertical conduit. Several specific forms are generally observed as defined in a system of flow regimes (Rouhani and Sohal, 1983) in conjunction with the different heat transfer regions affecting the liquid flow.
As the fluid enters the heated conduit at subcooled conditions, heterogeneous bubble nucleation begins to occur within the small pits and cavities which are designated as nucleation sites on the heated surface. These nucleation sites are activated when the temperature of the surface exceeds the saturation liquid temperature at the local pressure. The fluid flow remains predominantly single-phase. At a point depicted by the onset of nucleate boiling (ONB), bubbles begin to form and remain attached to the heated surface. As the bulk temperature of the liquid increases further downstream, the attached bubbles begin to grow and detach from the heated surface. Since the gas volume fraction remains low, the flow of individual ascending gas bubbles co-flowing with the liquid is referred as the bubbly flow regime. As the void fraction increases sharply at a location denoted as net vapor generation (NVG)-transitional between low void fraction region and region of which the void fraction increases significantly thereafter-a flow transition is exhibited whereby slugs of highly aerated liquid Fig. 1 Evolution of the steam-water flow in a vertical heated conduit (Yeoh and Tu, 2010; reproduced with permission © Elserier Ltd. 2010). move upwards along the conduit in the saturated nucleate boiling region. These so-called slug bubbles have characteristics of spherical cap nose and are somewhat abruptly terminated at the bottom edge. Slug pattern should be avoided since such two-phase flow structure causes undesirable flow instability. As large unsteady occurrence of gas volumes accumulates within these mixing motions, a flow regime known as churn-turbulent flow persists and the liquid may be flowing in an oscillatory fashion. At very high gas velocities, an annular pattern is observed whereby parts of the liquid flows along the pipe and other parts as droplets entrained in the gas flow. At the gas-liquid interface for sufficiently high gas velocity, there may be large amplitude waves that break-up during the flow process. The breaking of these waves is the continuous source of the deposition of droplets in the gas core. At even higher gas velocities, a disperse pattern exists in the form of mist and the liquid becomes severely deficient. There is now a considerable amount of liquid in the gas core and insufficient coolant in removing the heat produced by the nuclear reactor could lead to the melting of the core fuel element conduits.
In relation to PWR and BWR, it is imperative to note and recognize that the importance of the high-pressure range of these types of nuclear reactors. The pressure for PWR operations and transients is in the neighbourhood of 150 bar while BWR operations and transients is in the range of 50-70 bar. Liquid-to-gas density ratios at these two conditions are 6 and 18-27 when compared to a value of 1600 at 1 bar. These order of magnitude changes have important consequences not only on the dynamics of gas bubble growth (Patel and Theofanous, 1976;Theofanous et al., 1978) but also on the dynamics of the two-phase flow regimes.

Two-phase flow instability
Flow instabilities in two-phase flow processes are to be avoided particularly in water-cooled and water-moderated nuclear reactor systems including steam generators. Sustained flow oscillations in such systems may result in forced mechanical vibration of components or system control problems. The prevalence of flow oscillations that exist within the twophase flow processes can significantly affect the local heat transfer characteristics and induce the Boiling Crisis.
The classification of flow instability phenomena for water-cooled and water-moderated nuclear reactor systems can be derived from the considerations of the flow being subjected to static or dynamic instability. Table 1 depicts the different classes of flow instability that could be experienced in nuclear reactor systems. It is noted that a secondary phenomenon comes into effect after the occurrence of the primary phenomenon. This refers to particular cases when the occurrence of the primary phenomenon is a necessary Thermal hydraulic considerations of nuclear reactor systems: Past, present and future challenges 5 condition for the occurrence of the secondary phenomenon. Also, a flow instability is compound when several elementary mechanisms interact in the process and cannot be studied separately.
Analysis of instabilities occurring in natural circulation boiling loops is also very important for the safety of watercooled and water-moderated nuclear reactor systems (Prasad et al., 2007). Particularly for BWR, natural circulation is an important mode of operation for removing shutdown decay heat during accidents. In comparison to flow instabilities in forced circulation modes, flow under natural circulation is inherently less stable due to relatively small hydraulic driving head. In natural circulation loop, the heating process in the heater section is the driving force for the flow. The heat supplied will generate buoyancy and the flow will be created in the loop such that the buoyancy is balanced by friction in steady state. If the heating power is further increased, the flow rate will increase. On further increasing the heating power, the flow velocity may be so large leading to the insufficient time for the fluid to be heated and subcooled liquid enters the riser section and the buoyancy force reduces. The flow is decelerated and even reversed. A self-sustained oscillation is thus created which may be chaotic if the inlet subcooling is sufficiently large at a given heating power. These flow oscillations exist in two-phase natural circulation loop. In addition, BWR are subjected to control system and coupled neutronic and thermal-hydraulics instabilities. The former is due to the malfunction of reactor hardware. Suitable control mechanisms are required to be provided to handle this type of instability. The latter, also known as reactivity instability, is due to the void and power feedback effects on neutron kinetics and thermal-hydraulics respectively. In coupled neutronic reactivity instability, there is a power feedback (fuel transfer function) in addition to the flow feedback.

Understanding and mitigating potential severe faults
One special consideration on the design of nuclear reactors is the demonstration of safety under the circumstances of certain postulated severe faults. One hypothetical fault is the Loss of Coolant Accidents (LOCA) which may involve a double ended break of one of the main reactor recirculation lines. The important consideration of LOCA in nuclear design is due to its potential to evolve into a core melt accident under liquid deficient condition. Flow Boiling Crisis would occur in an LOCA under rapidly varying flow conditions. The desire to analyse such complex situations realistically brings about the specific need for a fundamental understanding of the Boiling Crisis. Ever since the Three Mile Island-2 (TMI-2) accident that occurred in 1979, the small break LOCA issue has attracted much attention in the nuclear safety. This subsequently resulted in the creation of the ECCS Rule 10 CFR 50.46 and associated Appendix K (Nusret, 2007). Theofanous et al. (1994Theofanous et al. ( , 1995, Asmolov et al. (2003) and Anh and Kin (2003) have focused on gaining a better understanding of the core melt accident progression in the reactor pressure vessel lower head and other related phenomenological uncertainties as a result of the TMI-2 accident. Core melt accident progression involves an early phase up to the partial melting of core material and a late Table 1 Classification of flow instability (Bourne et al., 1973)  phase resulting in the significant melting of the core material. During the early stage of the severe fault, the accident involves core uncover, heat-up, and partial melting in the reactor core region. The late phase is, however, characterised by significant core melting, core material relocation, and redistribution of predominantly melted core materials into the core region and reactor vessel lower plenum. The governing phenomena of the late phase melt progression involve a porous debris bed, molten pool and cavity formation due to the presence of very high temperatures, multicomponent and multi-phase materials, melting and freezing processes, and variable geometrical configurations.
If core melting proceeds, molten core materials (or corium being composed of conglomerated mixture of oxide component such as UO 2 and ZrO 2 , metallic component such as U, Zr, Fe and stainless steel, etc.) can accumulate on the core support plate and would be eventually relocated into the lower plenum region. In the presence of the lower plenum water, some portions of the relocating molten core materials will become fragmented into small solid particles in the lower head and the rest will maintain its original liquid phase. The main mechanism of fragmentation of the jet of flowing debris is a hydrodynamic process (Burger et al., 1995) unless there is a rather violent interaction approaching a steam explosion and the thermal interaction is a secondary contributor of the corium breakup. Rapid heat transfer from the jet of debris to the lower plenum water accompanies the hydrodynamic fragmentation process and debris oxidation, resulting in steam and hydrogen production and an abrupt increase of the primary system pressure.
Understanding the triggering of steam explosion in nuclear reactor systems is of major importance (Fletcher, 1995). Triggering can be initiated when a propagating wave is developed and led to a rapid transfer of heat from the melt to the water resulting in rapid rise of the local heat transfer and pressure. As observed experimentally, triggering is associated with the local collapse of the vapour layer around a melt droplet, followed by rapid fragmentation of the droplet (Corradini et al., 1988;Corradini, 1991). For core melt-water interaction, vapour film collapse may occur because water is forced into contact with the melt. This may be caused by an applied pressure pulse, the bulk flow of water or local coolant entrapment. Firstly, the pressure pulse induces a particle velocity in the coolant, towards the melt, at the liquid-gas interface. If this motion is adequate in driving the water into contact with the melt, triggering proceeds. Secondly, the bulk flow of water past a droplet causes the vapour layer to be advected away from the melt, causing the film to collapse. Thirdly, coolant being entrapped within the melt or against the vessel wall by the melt and is superheated until its temperature rises to the homogeneous nucleation temperature, at which point it flashes into steam being in contact with the melt, and triggering occurs. Explosions which result from a known trigger is referred to as triggered explosions and those occurring because of some uncontrolled event are known as spontaneous explosions. Explosion being artificially trigger, such as a detonator, is externally triggered.
In any nuclear safety analysis, the main challenge in mitigating potential severe faults or accidents is the reliability of thermal hydraulic computer codes in predicting the consequences for severe accidents and the specification of appropriate accident management strategies. It is therefore natural to suggest the following questions on the use of currently available thermal hydraulic computer codes: (i) which phenomena make different code predictions; (ii) which models and model parameters employed in a relevant code reflect the real physical phenomena; and (iii) which code is most appropriate for the code user using more of these codes to use. Answering all the above questions clearly is not simple because of: (a) the magnitude of model uncertainties and the scenario variability depending on the initial and boundary conditions that are usually not known a priori, (b) there are inherent difficulties in validating the models due to limited experimental data, and (c) the application of the models that are adequate for accident management strategies remain limited.

Flow regime identification and thermal hydraulic analysis
Flow regimes have been traditionally defined according to experimental visualisations carried out by viewing the twophase flow through transparent channels. Majority of all the reported data have been obtained through some methods of flow regime identification which employ a variety of signals and are based on two principally different approaches:
Thermal hydraulic analyses of nuclear reactor core can be performed using the safety analysis codes to estimate different thermal hydraulic safety margins. The safety margins and the operating power limits of nuclear reactor core under different conditions of primary cooling system, which is the system pressure, coolant inlet temperature, coolant flow rate and thermal power and its distributions are considered as the key parameters for nuclear safety analysis. Determination of flow transients through use of these codes is based on the solution of separate flow equations for the liquid and gas, known as the two-fluid model. Separate time dependent equations for mass, momentum and energy balance are employed for the liquid and gas in each element of space along the flow path. This allows for the prediction where the liquid and gas are not in thermal equilibrium, such as in the case of saturated vapor in contact with subcooled liquid during condensation, or the case of superheated steam carrying droplets of water at saturation temperature. The use of separate continuity and conservation equations in these computer codes makes it necessary to also include some equations in predicting the rates of exchange of mass, momentum and energy between the two phases. Appropriate forms of such equations, which are known as constitutive relationships, depend on the flow regime in actual situations.
A significant portion of nuclear reactor thermal hydraulic analyses have recently been performed via the multiphase computational fluid dynamic (CFD) codes such as ANSYS CFX, ANSYS FLUENT, STARCCM+, OpenFOAM, particularly the use of Eulerian-Eulerian, averaged two-fluid formulation in addressing different types of boiling regimes experienced in water-cooled and water-moderated nuclear reactors, both in normal operation and during design-basis and beyond design-basis postulated accident transients. In these CFD codes, separate time dependent equations governing the conservation of mass, momentum and energy are also solved for the liquid and gas phases, analogous to the two-fluid model in safety analysis codes. Appropriate constitutive relationships are utilised to model the interfacial rates of exchange of mass, momentum and energy between the two phases which they appear as source terms in the governing equations. In view of the simplistic models adopted in safety analysis codes, CFD can nonetheless describe the flow phenomena in much greater detail due to the increased resolution in mapping the entire flow configurations, leading to the solution in addressing selected nuclear reactor thermal hydraulic safety issues (Bestion, 2014), including the prediction of CHF.

Current state-of-the-art predictive models and methods for nuclear reactor systems
According to the thermodynamic definition of a phase, a two-phase flow refers to the situation when two phases are present and move simultaneously within the flow system. These two phases that co-exist simultaneously in the flow often exhibit relative motion among the phases as well as heat transfer across the phase boundary. In nuclear reactor systems, the type of two-phase flow being encountered typically consists of a mixture of water and air or of water and steam. Because of the complexity of the microscopic motions and thermal characteristics of the individual discrete constituents, solutions to the micro-level evolutionary equations can be prohibitive due to the uncertainty of these constituents at any instant of time and space. For most practical purposes of thermal hydraulic analyses, exact prediction on the evolution of details within the two-phase flow system is normally not required. Rather, the gross features of the fluid flow and heat transfer are of more significant importance. Owing to the complexities of interfaces and resultant discontinuities in fluid properties as well as from physical scaling issues, it is thereby customary to apply some sort of averaging process to the conservation equations, which leads to derivation of the conservation equations in the interpenetrating media framework (Yeoh and Tu, 2017).
In two-phase flow modelling, the averaging process that restricts the predictions to macroscopic phenomena allows the feasibility of solving the two-phase flow through suitable numerical techniques and the ease of comparison with the experimental data. Two-fluid modelling then proceeds by averaging the local instantaneous conservation equations for mass, momentum and energy over each phase. Numerous averaging approaches have been proposed. The details concerning a judicious choice of averaging can be found in Vernier and Delhaye (1968), Yadigaroglu and Lahey (1976), Delhaye and Achard (1976), Panton (1968), Agee et al. (1978), Banerjee and Chan (1980), Drew (1983), Lahey and Drew (1988), Besnard and Harlow (1988), Joseph et al. (1990), Drew and Passman (1999), Kolev (2005) and Ishii and Hibiki (2011). In general, averaging may be performed in time, space, over an ensemble, or in some combination of these. While averaging allows the mathematical solution of the problem tractable, there is a need to recover the lost information regarding the local gradients between phases, which can be re-supplied in the form of semiempirical closure relationships, also known as constitutive relationships. The success of a two-fluid model in predicting the two-phase flow heavily depends on the quality of the closure relationships that can be obtained for the various interfacial rates of exchange of mass, momentum and energy.

Critical considerations of the conservation equations for two-phase flow in safety analysis and CFD codes
Water-cooled and water-moderated nuclear reactors encounter two-phase mixtures of gas and liquid, either in their primary or the secondary coolant circuits under normal conditions, or everywhere during postulated accidental depressurisation. Safety analysis codes are designed to analyse the transient flow of such mixtures with a six-equation two-fluid model, where the two phases of gas and liquid have their separate equations governing the conservation of mass, momentum and energy.
Based on the critical review performed on the US NRC TRACE code by Wulff (2011), some of the many important limitations on the applicability of TRACE which can also be applied for other safety analysis codes to carry out thermal hydraulic analyses for nuclear reactor systems include: (i) Inherently developed for one-dimensional (1D) application due to the averaging of the conservation equations, these types of codes generally do not generally have a three-dimensional (3D) simulation capability because of: a) The fluid shear is generally ignored but is an essential element to predict counter-current flows. b) The use of fictitious body forces and fictitious distributed mass and heat sources to replace the contact forces at the wall, mass injection and wall heat fluxed for both phases. c) The need to impose perfect mixing thereby resulting in artificial damping and disables the ability of these types of codes to track the propagation of thermal and kinematic disturbances. d) The heavy reliance on numerical diffusion than the physical dissipation by turbulence. Numerical diffusion could reduce the ability of these types of codes to predict the onset of instability for nuclear reactors entering large power oscillations. (ii) The two-fluid model being employed is not closed for the different prevailing flow regimes in nuclear reactors and the consequences of that are unknown. Mixture model correlations being utilised in these types of codes bring about the questionable predictions of the phasic velocity and temperature. (iii) Mass may not be conserved and the models for boundary conditions are physically incompatible with conditions of boiling (especially in BWRs) and of flashing during Large Break LOCAs. (iv) Momentum may not be conserved because of incorrect averaging and boundary conditions. This is particularly important at branch points for the prediction of momentum transfer between, and the coupling in hydraulic networks by impedance and by inertia and subsequently the entire reactor component interactions. (v) Energy conservation equation may have been incorrectly averaged. Heat transfer from fuel cladding is not properly considered thereby leading to the incorrect consideration of clad conduction to accurately yield fuel pellet and clad temperatures whether for steady or transient, singlephase or two-phase flows. (vi) Boundary conditions for momentum and energy balances are restricted to flow regimes with single-phase wall contact. Limited constitutive relationships exist for solidfluid exchange of momentum and heat in prevailing flow regimes. For CFD codes, there are currently four different approaches that could be applied for the prediction of two-phase flows (see Table 2): (1) Averaging method, for example, the use of Eulerian-Eulerian two-fluid models which are analogous to the approach in safety analysis codes but with a 3D simulation capability.
(2) Bubble tracking method which is based on solving the equation of motion in the Lagrangian framework for each bubble and constitutive equations for external forces acting on the bubble. A two-way method is realised via representative source terms to affect the flow within the local instantaneous field equations in the Eulerian framework.
(3) Interface tracking method based on local instantaneous field equations and the use of auxiliary equations for tracking deformable interfaces. (4) Microscopic method which simulates the translational and collision of pseudo molecules (lattice gas method) or a molecular number density function (lattice Boltzmann method). As can be seen from the above description of the different methodologies, the spatial resolution of each method increases from (1) to (4) whereas its applicability to practical engineering problems decreases due to the need of accommodating the complex moving boundaries. Albeit the applicability of averaging and bubble tracking methods for practical problems, the accuracy of such methods comes into question because of the use of empirical formulations to resolve the constitutive equations. Despite the rapid progress in computer hardware and performance, interface tracking and microscopic methods remain applicable only to flows with bubbles that are confined in small-scale configurations. However, the ability of such methods to provide very detail information of the flow characteristics around bubbles greatly assists in the development of more accurate constitutive equations for the averaging and bubble tracking methods. A plethora of issues continue to exist within the CFD framework for two-phase flows and must be resolved to ensure that each method, as described in Table 2 with its own advantages and shortcomings, is reliable and accurate for the purpose of carrying out thermal hydraulic analyses for nuclear reactor systems.
3.2 Understanding some distinctive characteristics of two-phase flow to develop better predictive capability In LWR, HWR or PWR technologies, large-diameter flow channels are extensively used to increase the mass, momentum and heat transport capability of the working fluid. Most of the constitutive relationships employed in the active safety analysis codes have been developed based on the experimental data taken in small-diameter flow channels. It is thus vital that the field code development reflects the state-of-the-art research results of large-diameter channel flows in these constitutive relationships and analyse the uncertainty and scalability of these constitutive relationships in detail. Since it is often difficult to attain experimental data in high-temperature and high-pressure full-scale large systems, the knowledge obtained from the large-diameter flow channels under low-temperature and low-pressure conditions as proposed by Shen et al. (2018) is an important database which can be specifically utilised to approximately demonstrate the flow behaviours under high-temperature and high-pressure conditions through the original designed scaling methods. This will allow the capability of safety analysis and CFD codes to be extended to reliably predict flow regime transition in large-diameter flow channels.
Some of the pertinent characteristics of two-phase flows in large-diameter flow channels can be summarised in the following: (1) Owing to the interfacial instability of large Taylor cap bubbles, no stable slug flow persists in large-diameter channels. The flow regime transition from bubbly flow to cap/slug-type intermittent flow occurs gradually in large-diameter channels rather than the abrupt transition that prevails in small-diameter channels. Typical cap-turbulent flow in the large-diameter channel and its corresponding slug flow in the small-diameter channel are respectively illustrated in Fig. 2. A rather large liquid space around large cap bubbles in the large-diameter channels enables free movement of small bubbles, whilst small bubbles in the small-diameter channels are restricted by the channel inner wall and the large cap (or Taylor) bubble and these small bubbles are thus easily squeezed into the wake region of the proceeding large cap (or Taylor) bubble. This wake entrainment causes the slug (or Taylor) bubble to grow in length and greatly reduces the number density of the free small bubbles in the liquid phase in the small-diameter channel. The existence of the large cap bubble is usually confirmed by the observed local intermittent motions of small bubbles which move upwardly first then reverse violently. The whole shape of the large cap bubble cannot be identified clearly while it is travelling as visualised in the large-diameter channel.  (2) In both small and large diameter channels, wall-peaked void fraction profiles tend to occur at low void fraction and high superficial liquid velocities and core-peaked void fraction profiles tend to occur at high void fraction or low superficial liquid velocities. Both peaking heights are nonetheless significantly reduced in large-diameter channels, relative to those of small-diameter channels. This characteristic is illustrated by the comparison between the experimental data from small-diameter channels (Hibiki and Ishii, 1999;Hibiki et al., 2001) and large-diameter channels (Shen et al., 2006;Shen et al., 2018). Lateral movement of bubbles driven by intensive turbulence, thus significant secondary flow, in largediameter channels may account for this phenomenon (Ohnuki and Akimoto, 2000). (3) The occurrence of large cap bubbles rather than stable slug bubbles causes the differences between flows in large-diameter and small-diameter channels. It has been identified that traditional flow regime transition criteria such as Mishima and Ishii (1984) and others may not be applicable for large-diameter channels because of the observed flow regimes being different from the flow regimes in small-diameter channels. Hibiki and Ishii (2000) found that the flow regime transition from the bubbly to cap bubbly flows depends on the axial locations in the flow direction and the growing large cap bubbles increase the drift velocity and affect the void fraction in the vertical large-diameter channels. Schlegel et al. (2009) divided the flow regimes in vertical large-diameter channels into bubbly flow, cap-turbulent flow, churnturbulent flow and annular/mist flow. Their transition criteria are validated by air-water flow regimes identified by neural network and the flow regime data of Smith (2002) and Ohnuki and Akimoto (2000).
In the event of a severe fault such as an LOCA occurring in the primary circuit of the PWR, the depressurisation of the nuclear reactor results in steam generation. A favourable event that leads to the so-called reflux condensation causes the generated steam to flow into the steam generator through the hot leg. As steam condenses in the steam generator, the condensate will flow back through the hot leg to the reactor thereby resulting in counter-current steam/water flow. The success of the core being cooled depends to a certain extent on the behaviour of this counter-current flow. Nevertheless, the stratified counter-current flow of steam and condensate is only stable for certain ranges of steam and water mass flow rates. If the steam mass flow rate increases to a certain limit for a given condensate flow rate, a portion of the condensate may exhibit a partial flow reversal and be entrained by the steam in the opposite flow direction towards the steam generator. This phenomenon is known as the onset of flooding or the counter-current flow limitation (CCFL). Figure 3 illustrates the counter-current flow in the hot leg under reflux condensation conditions. For an additional increase of the steam flow, the condensate will be completely blocked and the reflux cooling mode ends. The cooling of the reactor core from the hot leg is impossible but may be continued by the coolant being drained through the cold leg to the downcomer. Figure 4 depicts the terminologies adopted to describe the counter-current gas-liquid two-phase flow in a model of PWR hot leg. The liquid film flows counter currently with the gas phase in hot leg channel for small gas flow rate. As the pressure difference remains low, and increases marginally with the air mass flow rate, this regime is denoted as the stable counter-current flow. With the gradual increase of the gas flow rate (m G ) a maximum gas flow rate persists in which the down-flowing water mass flow rate (m L,D ) in the reactor pressure vessel becomes equivalent to the inlet water mass flow rate. This point is defined as the CCFL. With the further increase of the air mass flow rate, the down-flowing water mass flow rate (m L,D ) approaches zero; this point is known as the zero liquid penetration (ZP). The region between the CCFL and ZP is defined as the partial delivery region. When the gas flow rate is decreased, a point is reached where a fully counter-current gas-liquid two-phase flow is established, which is the de-flooding point. Based on an extensive experimental database accumulated, Deendarlianto et al. (2012) comprehensively reviewed the development of phenomenological correlations and scaling parameters of the CCFL on the counter-current flow in a model of PWR hot leg. Most of the proposed correlations can nonetheless be applied under a relatively narrow range of conditions, generally limited to the test section conditions and/or geometry.

Review of models and methods for two-phase flow
One of the most practical models for the prediction of void fraction in gas-liquid two-phase flow is the use of drift-flux model (Hibiki and Ishii, 2003a). This widely-accepted model can be used to predict the interfacial drag as a constitutive model for the two-fluid model (Shen and HIbiki, 2013;Schlegel et al., 2017). However, the accuracy of the model predictions heavily depends on reliable correlations for the distribution parameter and drift-velocity for small-diameter as well as large-diameter flow channels. The drift-flux parameters are susceptible to inlet conditions for the cases of uniformly and non-uniformly introduced bubbles. Much work has been performed to develop the drift-flux correlations developed for small-diameter channels (Kataoka and Ishii, 1987;Hibiki and Ishii, 2003b;Shen et al., 2010a). To predict the flow behaviours in the transition from bubbly to cap-turbulent flow in long large-diameter channels, new drift-flux correlations have been developed by Shen et al. (2010b). Schlegel et al. (2013) overcome the weakness of the churn-turbulent flow prediction by improving the distribution parameter in this flow regime using the existing drift-flux correlations. Since there is no experimental data on annular flows in large diameter-channels are available due to the difficulty in achieving annular/mist flow in experiments, there is therefore no suitable drift-flux correlations that can be recommended.
In the quest of developing more physical-mechanism based models to predict the different flow regimes and transitions, transport equations for the Interfacial Area Concentration (IAC) and Average Bubble Number Density have been developed. This approach entails the component mechanistic modelling of various fluid particle interaction processes and is primarily developed to dynamically predict the interfacial transfer and the interfacial structure evolutions (Kocamustafaogullari and Ishii, 1995;Wu et al., 1998;Hibiki and Ishii, 2002;Yao and Morel, 2004;Cheung et al., 2007aCheung et al., , 2007bShen and Hibiki, 2013). To account for various 3D characteristic transport mechanism of cap-turbulent and churn-turbulent flow, the two-group transport equations for the IAC have been proposed (Smith et al., 2012). Bubble interactions with neighbouring bubbles and eddies are introduced in the transport equations through source and sink terms to predict the evolution of the IAC.
On a more sophisticated physical-mechanism based model, the homogeneous MUltiple-SIze-Group (MUSIG) model, which was first introduced by Lo (1996), has become widely adopted for predicting gas-liquid two-phase flow. In this model, the continuous bubble size distribution is approximated by discrete size fractions; the mass conversation of each size fractions is balanced by the inter-fraction mass transfer due to the mechanisms of bubble-bubble phenomena due to coalescence and break-up processes. The overall bubble size distribution evolution can then be explicitly resolved via source terms within the transport equations. Application of the homogeneous MUSIG model can be found in Pochorecki et al. (2001), Olmos et al. (2001), Frank et al. (2004), Tu (2004, 2005) and Cheung et al. (2007bCheung et al. ( , 2008. The inhomogeneous MUSIG model developed by Krepper et al. (2005) consisted of further sub-dividing the bubble phase into a finite number of velocity fields. This extension represents a robust and practical feature for two-phase flow modelling, particularly for bubbly flow simulations where bubbles can deform into different shapes and sizes. Useful information on the implementation and application of the inhomogeneous MUSIG model can be found in Shi et al. (2004) and Krepper et al. (2007a). In addition to the MUSIG model, Morel and Lavieville (2009) have applied a transport equation based on the conservation of the density S  of the moments of the bubble size distribution, which assumes to obey the log-normal probability distribution and first introduced by Lo and Zhang (2009) and more recently by Thakrar et al. (2015), to simulate boiling flows in a vertical pipe and in a vertical rectangular channel respectively.
With the aim of predicting the boiling process, the consideration of wall boiling models is crucial to the prediction of the bubble size distribution in the bulk two-phase flow. For two-fluid model, the Rensselaer Polytechnic Institute (RPI) boiling model from Kurul and Podowski (1990) represents the cornerstone for the prediction of boiling flows through its use in safety analysis and CFD codes (Tu and Yeoh, 2002;Tu et al., 2005;Tu, 2006a, 2006b;Yeoh et al., 2008;Krepper et al., 2013;Colombo and Fairweather, 2016). Principally, the RPI boiling model considers the heat flux from the wall being partitioned between the mechanisms responsible for the heat transfer process: single-phase convection, quenching and evaporation. Even the RPI model was conceived to be built mechanistically, almost all the RPIbased models have been forced to rely on some empirical or semi-empirical closure relation, particularly for the evaporative heat transfer contribution, which requires knowledge of the active nucleation site density, bubble departure diameter and frequency. In general, poor predictive accuracy of these models has been achieved when compared against subcooled boiling data over a wide range of mass and wall heat fluxes, and inlet subcooling. A comprehensive evaluation of existing correlations has been performed in Cheung et al. (2014) to assess the performance of these empirical models. To circumvent the deficiencies with the RPI-based models, the wall heat flux partitioning model was improved by considering the likelihood of sliding bubbles at the heated wall yielding an additional heat flux, re-determine the bubble departure diameter and frequency via the force balance model and fundamental consideration of the growth and waiting time and the fractal approach for the active nucleation site density. More details on the salient features of the model can be found in Yeoh et al. (2014).
In any nuclear reactor systems, thermal hydraulic analysis of reactor core is normally performed via the use of subchannel analysis code to estimate the various thermal hydraulic safety parameters in the event of an FBC such as Critical Heat Flux (CHF) ratio, Critical Power Ratio (CPR), fuel centre line temperature, fuel surface temperature, subchannel maximum temperature and bulk coolant outlet temperature (Chelemer et al., 1977). CHF ratio and consequently, the CPR and fuel centre line temperature are the main parameters limiting the maximum operating power of the reactor (Cheng and Muller, 2003). Reactor core fuel assemblies are of different shapes: square, hexagonal or circular in cross section such as illustrated in Fig. 5 (Ginox, 1978;Todreas and Kazimi, 2001). Typical PWR and BWR fuel pins are assembled in the form of regular arrangements/ patterns and are known as fuel assemblies/bundles. Inside the fuel assembly, fuel pins can be arranged either in square or triangular lattice configuration.
Thermal hydraulic safety analyses of nuclear reactor are performed in two ways. Firstly, the so-called system level safety codes such as RELAP, TRAC, TRACE, RETRAN, CATHARE, CANAC, ATHLET, NOTRUMP, etc. are utilised to obtain the system behaviours under different steady state and transient operating conditions. Detailed analysis of the reactor core is nonetheless performed using the sub-channel thermal hydraulic codes such as HECTIC, ENERGY, SUPERENERGY, COBRA (-I, II, IIIC, IV), CANAL, HAMBO, FLICA, THINC, VIPRE and COBRA-FLX (a comprehensive review of these sun-channel analysis codes can be found in Moorthi et al. (2018)). A sub-channel can be defined as a flow passage formed between number of rods or some rods and wall of channel/shroud tube. The sub-channels can be formed by either coolant centred sub channels or the rod centred sub-channels as shown in Fig. 5. The concept of sub-channel analysis method is an important tool for predicting the thermal hydraulic performance of the local conditions of rod bundle nuclear fuel element under single-phase as well as two-phase. The analysis considers a rod bundle to be a continuously interconnected set of parallel flow sub-channels that are assumed to contain one dimensional flow coupled to each other by cross flow mixing. The typical approach for sub-channel analysis with increasing complexity is depicted in Fig. 6.
On the use of CFD codes for sub-channel analysis, CFD simulations have been performed for PWR fuel assemblies to estimate the single-phase velocity and temperature fields , secondary flow in rod bundles , turbulent flow structures and mixing across the narrow gap between rod bundle sub-channels (Chang and Tavoularis, 2008). Various CFD models have been adopted including RANS (Chandra et al., 2009) or URANS (Merzari et al., , 2008 with different turbulent models k-ε (Baglietto and Ninokata, 2005), LES (Biemüller et al., 1996) and DNS (Ninokata et al.,  Thermal hydraulic considerations of nuclear reactor systems: Past, present and future challenges 13 2004; Baglietto et al., 2006). Applications of CFD for the prediction of two-phase flow velocity, sub-cooled boiling for void distribution in Krepper et al. (2007b). Analyses extended in the rod bundle (Carver et al., 1984;Anglart and Nylund, 1996). The CFD approach permits the feasibility of modelling complicated flow field consisting of a spacer grid and a rod bundle and to evaluate the local velocity and enthalpy distribution around the rod surface, which can be taken as the initial conditions for the onset development of the two-phase structure or for improved CNF performance.
To better understand the dynamics around the CCFL in a PWR hot leg, several types of analytical studies have been realised, mainly focussing on the prediction of the gas velocity at the inception of flooding, and the assessment of the proposed analytical models to the real size of PWR hot leg. Available analytical studies are: (a) 1D stratified twofluid model, (b) semi-analytical and (c) scaling parameters.
For 1D stratified two-fluid model, three different concepts have been proposed to determine the CCFL in a PWR hot leg. Based on the observed flooding mechanisms of the experimental work of Siddiqui et al. (1986), Ardron and Banerjee (1986) developed a model to solve the steady state mass and momentum balance equations for a countercurrent stratified flow of a liquid film and a gas in a horizontal pipe. For this first concept, the CCFL is assumed due to the onset of slugging in the lower leg of the elbow close to the bend, where the liquid depth is greatest. The liquid level from the hydraulic jump to the horizontal leg outlet decreases continually. A free out-fall is assumed at the horizontal leg exit, implying that the liquid velocity at this point is equivalent to the critical velocity. The envelope theory which represents the second concept is applied to solve the momentum equations in the two-fluid model. By assuming uniform flow in steady state, no entrainment, and no phase change, the envelope is employed to the locus of tangents to the operating lines in the superficial gas and liquid velocities plane for a constant void fraction and therefore represents a limiting curve separating the operating region from an unattainable region for counter-current flow (Bankoff and Lee, 1985). Based on this concept, CCFL is considered to correspond to the retardation of the liquid film due to the interfacial shear stress. Models based on this concept include: Ohnuki et al. (1988), Geweke et al. (1992), Minami et al. (2008), and Nariai et al. (2010). The third concept implements the idea on the onset of slug flow criteria as a flooding criterion in a model of the PWR hot leg. Here, the flooding point is determined by the integration of the momentum equations of a 1D two-fluid model, whereas the slug criterion has been employed as a boundary value. Introduced by De Bertodano (1994), the important parameters in this concept are three friction factors and momentum loss of the elbow. Hydraulic losses at the elbow are modelled by extending the length in the horizontal section. Kawaji et al. (1991) developed two semi-analytical methods to determine the CCFL in the complex piping system. Slugging mechanism near the elbow and entrainment of liquid droplet has been considered. In order to accurately determine the relative velocity between the gas and liquid near the elbow at the flooding inception, it has been assumed that the annular flow in the vertical pipe is maintained for some distance past the elbow in the entrance section of the inclined leg. The gas void fraction near the elbow can subsequently be calculated from the available liquid film thickness equation of annular flow. For the first mechanism being the slugging mechanism near the elbow, Kawaji et al. (1991) modified a correlation of slugging criterion near horizontal pipe proposed by Taitel and Dukler (1976). For the second mechanism, Kawaji et al. (1991) considered the entrainment and carry-over of droplets generated as a result of the breakup of the turbulent jet-like liquid stream. The consideration of break-up is necessary but not a sufficient condition for flooding as the droplet must be entrained and carried upstream by the gas flow. A simple force balance on a droplet of given diameter moving in a counter-current flow of a gas with a relative velocity was considered. In the case of analytical or semi-analytical works on the CCFL in a PWR hot leg using multiple elbows was found to be limited. For this purpose, Kawaji et al. (1993) adopted a superposition principle whereby the piping system is represented as a combination of simple geometries (vertical, horizontal or inclined pipes). The most limiting gas velocities as each liquid flow rate are combined to obtain the predicted CCFL conditions over the entire range of the liquid flow rates. Glaeser (1992) proposed the determination of suitable parameters to correlate the CCFL data from the view point of analytical work as a possible solution. A simple theoretical basis for the extension of the flooding equation was provided in order to decide the suitable scaling parameters between the Wallis parameter and the Kutateladze number of counter-current flow in a full geometry reactor scale. In this work, the scaling parameter in a separated counter-current flow in a PWR hot leg is assumed for the configuration of a horizontal pipe.
In recent times, detailed 3D information of the transient behaviour around CCFL in a model of PWR hot leg has become a new requirement in reactor safety analysis. CFD allows substituting geometry-dependent empirical closure relations with more physically justified closure laws that are formulated at the scale of the structures of the gasliquid interface. Introducing CFD to resolve the countercurrent flow in a model of PWR hot leg includes the investigation of CCFL mechanisms, heat transfer effects, flow patterns, hysteresis behaviour, and the extension of the obtained flow behaviour from small scale to full reactor scale. The serial works of Wang and Mayinger (1995), Murase et al. (2009 and  and Kinoshita et al. (2009) are noted. Here, we review the development of a general model being closer to physics and including less empiricism by the Deendarlianto and his co-workers in Dresden, Germany.
As described in Höhne (2009) The above result  indicates that the AIAD model is a promising way to simulate the phenomena around CCFL in a PWR hot leg. Moreover, further improvements are required. The usage of the morphology detection algorithm should also be made possible to resolve vertical flow regimes. Therefore, it would be necessary to include the modelling of non-drag forces (lift force, wall lubrication force, virtual mass force, etc.) in the AIAD model as well as the available complete for polydispersed flows. Turbulence damping procedures should include the existence of small surface instabilities in the macroscopic model. In the last decade, the increasing use of 3D CFD codes to explain the phenomena around CCFL has surpassed the use of 1D system codes with the required accuracy and spatial resolution. CFD codes contain models for simulating turbulence, heat transfer, multiphase flows, and chemical reactions. Nevertheless, validations are still required to be performed for the these CFD codes to achieve a high level of confidence to carry out thermal hydraulic analyses for nuclear reactor systems. Attainment of more experimental data to assess the CFD models remains paramount in order to obtain more reliable models.

Review of measuring methods for two-phase flow
Validation of predictive models for two-phase flow requires an accurate measurement of the void fraction and inter-facial area concentration in nuclear reactor systems. A variety of instrumentation is available for measuring the void fraction. In general, void fraction can be measured using capacitance, conductance, optical, ultrasonic, or radiative methods (Boyer et al., 2002).
Radiative methods such as gamma ray attenuation and high energy X-ray methods are non-intrusive chordal, xy-tomographic, or yz-tomographic measurements of the void fraction; z is taken to be the flow direction. If the void fraction is a strong function of time or geometrically skewed in the pipe (e.g., flow regime effects), chordal measurements can generate more uncertainty (Abro and Johansen, 1999). Multiple detectors are required to obtain accurate measurements in the case of varying flow regimes. XY-tomography measures through the flow plane and can provide an Eulerian view of the flow as it passes through the detection region via post processing of the data with techniques similar to those used for CT and MRI scans in the medical fields. Resolution of these images is dependent on the number of available detection angles to map the two-phase flow. Ultrafast X-ray systems, which utilise an X-ray beam that rotates about a static target along a wave guide with a co-located detection ring, can allow for higher resolution and scanning frequency in a smaller physical space than typical radiative measurement systems (Banowski et al., 2016). The second tomographic imaging technique takes a two-dimensional (2D) density picture along the flow direction and one transverse direction. This method suffers, however, from geometric problems such as curved pipe surfaces that could lead to distortion of the image and assumes that the image is attained at a sufficiently fast resolution with the flow being static during the imaging (Banowski et al., 2016).
Another non-intrusive technique-ultrasonic tomographyemploys sound instead of radiative means in measuring the void fraction. Emitted ultrasonic pulses are recorded with various receivers located within the test section. Large number of sensors are required to be positioned around the piping. In order to obtain an accurate representation of the void fraction. This method has limited applicability for low gas fractions due to a loss of the ultrasound transmission, thereby providing incorrect results, and require significant amounts of post-processing to view images, so on-line viewing is restricted to approximately 10 Hz. Data can nonetheless still be collected at much higher rates and postprocessed for void fraction measurements. The processing of ultrasonic tomography requires significant computational effort and small changes in the setup geometry can have a large impact on the accuracy of the results (Rahim et al., 2007).
The use of optical probes to measure void fraction is restricted to only point measurements. These probes generally operate based on Snell's reflection law for optics. As different phases pass over the probe tip the index of refraction changes, this results in a change in the voltage measurement from a photodiode. Since the probe is a single point measurement like a thermocouple, data is analysed statistically over a period of time to generate the timeaveraged multiphase parameters (e.g., void fraction, interfacial area concentration, and interfacial velocity). Being the nature of the device, optical probes are unable to determine bulk parameters unless operating in steady state conditions and multiple measurements are required to be strategically placed at different locations (Chabot et al., 1998;Choi and Lee, 1990).
Electrical methods, which rely on the difference in the permittivity or conductivity between the fluid phases to determine the presence of voids within the two-phase flow, can exist in the form of a single probe like the optical probes, a large conductance cell across the pipe, or an intrusive method capable of direct tomographic imaging by introducing a grid of sensing points that allow for full reconstruction of the area void fraction. Point probes suffer from the same issues as aforementioned. Individual probes have been developed with up to four measuring points in order to evaluate bubble size, contact angle, and velocity (Smith, 2002;Sun et al., 2002;Shen et al., 2006Shen et al., , 2018Schlegel et al., 2012). For 3D two-phase flow measurement, the four-sensor probe consists of one central front sensor (0) and three peripheral rear sensors (1, 2, 3), which are usually made from either optical fiber (optical type) or metal needle (conductivity type) as shown in Fig. 7. These sensors detect the time when the phase changes from one to the other. Then the conversion of the known distances between sensors and the detected time into local flow parameters can be reached through mathematical processing algorithms. These algorithms enable the four-sensor probe to indirectly measure the important local flow information. However, their accuracy is dependent on the algorithms used to relate the measured value to void fraction, which can change significantly with flow regime (Lee et al., 2017). The grid method, formally known as wire-mesh sensors (WMS), measures the void fraction in multiple small conductance cells distributed across the pipe cross-section, but also perturbs the flow that is not present with radiative and ultrasonic methods. WMS directly visualise the flow for measurement of the void fraction as well as for measurements of bubble velocity and size. Measurement of phase with WMS is carried out by placing two perpendicularly oriented planes of parallel electrodes separated by a small gap in the direction of the flow (see Fig. 8). One plane will act as a transmitter, while the other will be the receiver. The transmitter plane yields a 6 micro-second bi-polar voltage pulse on each electrode. A potential is built up between the driven transmitter and the receiver wires, which causes a current proportional on the resistivity of the fluid between the electrodes, to flow to the receiver electrodes. This change in current is subsequently measured by a transimpedance amplifier and ultimately passed to an analog-to-digital converter (ADC). The magnitude of the pulse being measured is proportional to the conductance of the mixed phase between the electrodes. Each electrode in the transmitter is pulsed individually to isolate the current source location and a tomographic image of the flow is obtained as seen in Fig. 9.   Fig. 8 Schematic of wire-mesh sensors. Transmitters on the left supply a current to the electrodes sequentially and are measured by the receivers on the right (Prasser et al., 1998; reproduced with permission © Elsevier Science Ltd. 1998).

Fig. 9
Wire mesh sensor timing for voltage pulses during operation and subsequent sampling. Labels are defined in the previous Fig. 8 (Prasser et al., 1998; reproduced with permission © Elsevier Science Ltd. 1998).

Future modelling challenges in the design of advanced nuclear reactor systems
For the future of thermal hydraulic applications, CHF remains the most important parameter in the design, operation and setting safety limits of nuclear reactors for maximum efficiency and maximum power throughput. Understandably, the phenomena and details of the CHF mechanisms are very complex and diverse. Whether for designs of current and advanced nuclear reactor systems, a full description of CHF (under all possible situations) still presents many challenges especially in the reliable thermal hydraulic prediction of CHF (Yadrigaroglu, 2014). Other future modelling challenges in the thermal hydraulic considerations of nuclear reactor systems include understanding the potential use of nanofluids for high boiling heat transfer and CHF enhancement, the internal flow induced vibration mechanism in triggering the FBC, the ocean motions of small modular reactors on floating platforms and the singlephase and two-phase flow and heat transfer of other types of coolants in Generation IV reactor technologies as compared to the current water-based nuclear reactor systems.
4.1 Current capability of safety analysis, subchannel analysis and CFD codes to predict CHF for water-cooled and water-moderated nuclear reactor systems CHF is essentially the inability of the wall to evacuate the heat flux imposed, take for an example the case of a nuclear fuel rod in the nuclear reactor core resulting in an overheating condition of the wall such as the fuel rod cladding at high heat flux. Under normal flow conditions, there are basically two broad categories of CHF. During FBC, NVG will take place at low quality, under typically nucleate-boiling and bubbly-flow conditions by some overcrowding of bubbles or starvation from liquid of the heated wall. In the case of dryout which is expected under high-quality, annular-flow conditions, depletion or drying out of the liquid film on the heated wall produces the CHF condition (Tong and Hewitt, 1972). The industrial challenge regarding CHF in nuclear reactor systems is the ability to accurately predict the CHF condition in LWR fuel bundles. CHF governs the performance of the fuel and of the reactor core under normal operating conditions. Some or all the difficulties of the classical CHF prediction methods for fuel bundles include the applications of empirical correlations and/or subchannel analysis and the increasing use of CFD to simulate the flow in a reactor fuel element bundle and simultaneous predict the CHF condition in a more mechanistic manner.
Large heat transfer coefficients exist during subcooled flow boiling. This enhanced heat transfer mechanism is nonetheless limited by the value of the CHF. The coefficient of heat transfer decreases in the case of dryout and this leads to a rapid excursion of the heater temperature which results in heater destruction due to melting. Many empirical correlations for CHF have been obtained, and they are currently applied in purpose-developed 1D safety analysis codes. These correlations are nonetheless limited to a specified region of fluid conditions, fluid properties, and specific defined geometry (Krepper et al., 2007b). An extension to the fluid parameter validity was overcome by use of lookup tables based on experimental data. The vast amounts of available CHF data (tens of thousands of points) were collected and assembled in data bank files; ways of interpolating the data to find the desired result for the specific condition of interest were established (Doroshchuk et al., 1975;Groeneveld et al., 1996Groeneveld et al., , 2007. These system analysis codes generally adopt the quasi-steady-state assumption for all closure laws (i.e., the correlations are evaluated based on the instantaneous values of their parameters, ignoring any transient effects) and the prediction of CHF. Although steady-state CHF has been studied extensively, much less is known about the occurrence of CHF under transient conditions, such as LOCA conditions or during transients. The difficulty on how to predict these under transient conditions is akin to the of prediction of heat transfer from wall to coolant under transient flow conditions. In both cases, all the correlations and other closure laws upon which the solutions are principally applied only to steady state condition. The outcomes certainly depend on the very nature of the transient and cannot be generically described, and even less, correlated. A couple of rare analyses of the transient CHF problem can be found in Pasamehmetoglu and Gunnerson (1985) and Pasamehmetoglu et al. (1987).
Subchannel analysis can be performed to: (a) predict accurately the flow conditions in different parts of the bundle cross-section (the subchannels) using a subchannel code which accounts automatically the heat flux distribution, and (b) link and correlate local flow conditions to a local CHF criterion. This arises nonetheless the need for a subchannel-conditional CHF criterion or model that is not easily to be obtained directly. Use of lookup tables for a simple tube configuration or extracting such local CHF correlations from experimental data remains the only feasible way of carrying out the subchannel analysis. Accuracy of the subchannel analysis depends on the adequacy of the closure laws used to describe the wall and interfacial relationships and the transfers between subchannels, which constitute the many challenging issues being encountered in the current state-of-the-art analysis employed in the nuclear industry (Gluck, 2007). Subchannel analysis also shares the same difficulties with the 1D system analysis codes on the prediction of CHF under transient conditions.
The use of CFD codes by their own nature focuses on the basic considerations of turbulence and transient conduction but these codes still adopt empirical closure laws, for example, in describing the forces acting on bubbles. CFD codes are inherently not limited to steady state condition and should be more applicable to treat difficult transient situations in nuclear reactor systems. All the difficulties confronted by the classical CHF methodscorrelations and/or subchannel analysis-would have been eliminated by CFD simulations of the flow in a bundle and simultaneous prediction of the CHF condition mechanistically. The phenomena being prevalent in nuclear reactor systems are multi-scale in nature. At the microscopic scale of the bubble, the activation of a single and of a collection of nucleation sites, the growth, detachment and possible sliding of individual bubbles, bubble interactions, bubble coalescence and crowding near the hot surface, creation of a hot patch (conjugate heat transfer with the wall) and overheating of the wall, are required to be considered. At the scale of a bundle of subchannels, the distribution of voids and other flow variables become important. At the scale of the bundle, the exchanges between subchannels, including any effects of spacers, channel wall and other aspects of a fuel bundle should be determined to predict the flow behaviour along the bundle. Clearly the problems at each scale will have to be tackled with the appropriate tools to predict the phenomenon. These generally involves the use of: (a) interface tracking method to predict the growth of a single boiling bubble on the heated wall (Lay and Dhir, 1995;Son et al., 1999), coalescence of bubbles from multiple nucleation sites and blanketing of the heater surface with vapour (Son and Dhir, 2008), (b) two-fluid models and N-phase mixture model with algebraic slip between the phases and a temperature equation for each phase to predict the axial void fraction profiles near a heated wall (Thomas et al., 2013), (c) high resolution mesh in conjunction with the two-fluid model to resolve the flow and predict the conjugate heat transfer from the rods of a PWR 5  5 rod bundle (Dominguez-Ontiveros et al., 2012), and (d) Eulerian-Lagrangian approach to simulate steamdroplet flow in a heated annular-mist flow regime to investigate the phenomena of droplet deposition (Caraghiaur, 2012;Anglart and Caraghiaur, 2011) and to account for the complex geometrical effects of the spacer grids.

Heat transfer enhancement via nanofluids in nuclear reactor systems
Nanofluids, which consist of improving the heat transfer characteristics of fluids by the addition of solid particles with diameters below 100 nm, are being considered as potential working fluids to be used in high heat flux systems such as nuclear reactor systems. Current findings as reviewed by Wu and Zhao (2013) to the suitability of application of nanofluids in energy-efficient heat transfer systems are: (a) The flow and heat transfer performance of thermal systems greatly depend on the understanding of key thermo-physical properties of fluids. Many experimental investigations have nonetheless revealed that not all nanofluids with enhanced conductivity can bring about forced heat transfer enhancement. Under suitable conditions, a higher heat transfer coefficient of nanofluids can be achieved compared with the base fluid with the appropriate nanoparticle material and size, particle concentration range, and preparation method and additives to maintain stability. (b) In general, nanofluids behave like single-phase fluids conditionally, depending largely on the base fluid, nanoparticle materials, concentration and size. Under forced convective flow condition, either in the laminar or turbulent regime, the shear stress and temperature gradient in the boundary layer may re-distribute the nanoparticles, resulting in a non-uniform thermal conductivity and viscosity distribution thereby altering the thermal resistance of boundary layer. Migration or movement of nanoparticles from Brownian motion or thermophoresis should also be determined if the nanofluids behave like two-phase flows in which the slip velocity between the particle and base fluid plays an important role in heat transfer performance. (c) Pool boiling encompasses aspects of bubble generation, growth, and detachment from heated surfaces, motion and coalescence in the fluid, and bursting at the liquid surface. Nanoparticles suspended in nanofluids could affect the bubble growth and detachment from heated surface, motion and coalescence in fluid. The presence of nanoparticles should result in a complete change of flow field and temperature distribution of boiling nanofluids, and on the boiling heat transfer for the number of generated bubbles, bubble growth rates, the size and frequency of detached bubbles. State-of-the-art dynamic particle image velocimetry by Dominguez-Ontiveros et al. (2010)  and CHF enhancement of nanofluids can increase the peak-cladding-temperature margins (in the nominalpower core) or maintain them in uprated cores for ECC in PWRs during a large-break LOCA. Bringing the decay heat from the molten core out more efficiently can increase the margin for vessel breach by 40% for the IVR in advanced LWRs during core molten accident (Buongiorno et al., 2008;Pham et al., 2012). Given the preceding findings, there remain many challenges ahead to better understand the underlying mechanisms leading to the unprecedented thermal transport phenomena. Long-term physical and chemical stable nanofluids with non-agglomerated nanoparticles for highvolume production are required. Anomalously high transport properties of nanofluids need adequate resolution. Future works in heat transfer enhancement of nanofluids need to be focused on: (1) Development of nanofluid stability under both quiescent and flow conditions. Such a property determines whether nanofluids can be applied successfully in engineering systems where steady thermal management is required.
(2) Creation of a nanofluid database of thermo-physical properties, including detailed characterization of nanoparticle sizes, distribution, and additives or stabilizers. Priority should be directed towards nanofluids with promising potential. (3) Extensive experimental and numerical studies on the interaction of suspended nanoparticles and boundary layers should be performed to uncover the mechanism behind convective heat transfer enhancement by nanofluids. Numerical simulation methods such as Lattice Boltzmann dynamic simulations could provide the necessary insights into the motion and distribution of nanoparticles in the boundary layer of fluids. Promising nanofluids that could boost the convective heat transfer increase, without much concomitant viscous pressure drop increase, need to be also considered. (4) Bubble dynamics of boiling nanofluids should be investigated experimentally and numerically. By considering the surface tension as well as influences of the presence of nanoparticles and additives if used, exact contributions of solid surface modifications and suspended nanoparticles to CHF enhancement in boiling heat transfer could be identified. Lattice Boltzmann simulations to the numerical analysis of bubble motions as suggested by Takada et al. (2000) could provide vital information to CHF enhancement as such a methodology possesses the relative ease of boundary condition implementations on complex geometries, high efficiency on parallel processing, and flexible reproduction of interfaces between phases.

Flow-induced vibration in nuclear reactor systems
Two-phase Flow Induced Vibration (FIV) can exist in heat exchangers as well as in nuclear power plant components (Konno and Saito, 1985;Pettigrew and Taylor, 1994;Miwa et al., 2014). In avoiding structural damage that may be caused by fluid-solid interaction, knowledge of two-phase gas-liquid flow induced force fluctuation magnitudes and its predominant frequency is paramount. Separate models for the fluid and structural dynamics are normally developed to predict the FIV phenomena, and thereafter coupled through the hydrodynamic and structural force terms. In most cases, models available for structural dynamics are modelled as linear-oscillator. However, thermal hydraulic models for fluid dynamics are generally more complicated since the fluid flow is inherently nonlinear and possesses multi-degree-of-freedom behaviours. In addition, mechanisms of two-phase FIV can be quite different from single-phase flow, due to complex motions/interactions at phase boundary, differences in material properties (density, surface tension, viscosity, etc.) and phase change process via energy transfer/generation. There are essentially three categories for two-phase FIV, which are based on the JSME handbook (JSME, 2003). The first category being momentum fluctuation, is caused by the density difference between two phases, and large FIV is generated due to the change in flow direction and the impact force of two-phase mixture acted on piping component structure such as elbows and T-junctions. The second category being thermal-hydraulic vibration associated with phase change, is induced from the nature of two-phase flow which involves phase change due to the energy transfer through interfacial boundary and/or energy generation within the phase such as boiling and condensation, which will involve highly unstable and oscillatory behaviours, and FIV is promoted within two-phase flow systems. The third category being bubble-induced vibration, is due to the dynamics of various shapes and sizes of bubbles that induce sloshing, fluctuations and disturbances within the flow fields. All the FIV categories are coupled with the fluctuating characteristics of two-phase flow, namely momentum, pressure and void fraction fluctuations. These phenomena are depended on the flow orientations: axial flow, internal flow, and cross flow.
Vibration mechanisms such as fluid-elastic instability, hydraulic oscillations due to phase change, random turbulence excitation and acoustic noise are the primary vibration excitation mechanisms in axial two-phase FIV (Taylor and Pettigrew, 2001). Such applications vibration phenomena are prevalent in BWR fuel rods and tube bundle of the steam generators for PWR (Pettigrew and Taylor, 1994;1998, Feenstra et al., 2009Chu et al., 2011). The effect of void fraction fluctuations is the significant two-phase flow parameter for axial FIV phenomenon. Studies show that the damping effect becomes important for the void fraction less than 80%. For cross flow two-phase FIV, primary causes of the FIV phenomena are due to the pressure fluctuations generated by Karman vortex, collisions of the bubbles against structure surface, and turbulent eddies interacting against structures. Additional three key mechanisms cross flow two-phase FIV have also been identified by Pettigrew et al. (1998) which include fluid-elastic instability, periodic wake shedding and random excitation due to turbulence. The occurrence of these phenomena depends on the flow regime, which is associated with physical properties of fluids, inlet flow conditions, local void fraction, and flow channel geometries. Cross flow two-phase FIV exists in shell-and-tube heat exchangers (Mitra et al., 2009;Sasakawa et al., 2005;Zhang et al., 2007;Kanizawa et al., 2012) and U-tube region in PWR steam generators (Pettigrew et al., 1998). With regards to internal two-phase FIV, flow turning elements are considered as one of the major sources for causing FIV by generating sudden change in momentum flux, pressure fields, or creating secondary vortices due to boundary layer separation Cargnelutti et al., 2010;Pontaza and Menon, 2011;Yamano et al., 2011). The response of the structure typically increases linearly (in time) when the dominant fluctuating frequency overlaps with natural frequency of the structure. Such phenomenon is known as the resonance and should be avoided to ensure safe system operation.
Much progress has been accomplished to understand the FIV mechanisms particularly for external two-phase FIV such as axial and cross flow (Miwa et al., 2014). Nevertheless, the thermal hydraulic predictive capability for internal two-phase FIV is still under developed. The vibration of flow channel significantly influences the two-phase flow structure. From the two-phase flow perspective, majority of the observed phenomena are flow regime specific. In general, slug and churn flow regimes have the strongest fluctuation compared to other flow regimes. Effect of twophase flow regime on the force needs to be investigated, particularly for large diameter pipe two-phase flow and annular two-phase flow regime. Possible considerations of the geometric and operation conditions could minimise internal two-phase FIV especially by raising the natural frequency of the piping structure than the two-phase characteristic frequency and minimise the occurrence of liquid slug when operating under two-phase flow condition. Impact force of liquid slug against piping wall can usually be considerably large, resembling like the phenomena of water-hammer.

Nuclear reactor systems due to ocean motions
Recently, the development of advanced small modular reactor (SMR) is attaining significant traction. According to the classification adopted by the IAEA, SMRs are reactors with the equivalent electric power less than 300 MW (IAEA-TECDOC-1451. Compared with the Generation III nuclear reactor technologies, the SMR of modular design is more flexible and of higher passive safety. SMRs being of interest for seawater desalination, hydrogen production and coal liquefaction could thus be built on land or on a floating platform. A new concept for offshore nuclear power plant (ONPP) with enhanced safety features has been proposed by Lee et al. (2013). The design concept includes the mounting of the nuclear reactor systems on gravity-based structures and a new emergency passive containment cooling system and emergency passive reactorvessel cooling system. The ONPP combines the state-ofthe-art LWR and floating platforms similarly used in offshore oil and gas operation. Since seawater was used as coolant against earthquakes, tsunamis, storms, and marine collisions, its ocean-based passive safety systems eliminate the loss of ultimate heat sink accident by design.
Design and operation of nuclear reactor systems floating on the ocean is different from that on land. Owing to the action of waves, the coolant flow in primary loop and secondary loop will be affected by additional forces. The gravitational pressure drop may also be altered due to the height variation between the reactor core and steam generator. The ocean motions contributing to the nuclear reactor thermal hydraulic characteristics are caused by the ocean surface waves that occur in the upper layer of the ocean. They usually result from wind or geologic effects and range in size from small ripples to huge waves and thereafter affect the ship or floating platform motion. The law of ocean waves is complex. Wind usually changes with space and time. Therefore, the ocean motions and ship or floating platform motions change with space and time, just like the wind waves. The ship or floating platform motions can include heeling, heaving, rolling, pitching, yawing, swaying and surging motions (Ishida et al., 1990). In heeling or inclination motion, there is only height difference change but no additional force, which is different form the other six motions. The amplitudes of yawing, surging and swaying motions are usually very small, since the ships and floating platforms are long and narrow designed. Compared with these three motions, the heeling, heaving, rolling and pitching motions are the most prevalent motions (Yan, 2017).
In heeling motion, neither additional acceleration nor inertial force exists. Only the relative height difference between reactor core and steam generator is changed. This represents the simplest motion and closest to the stationary state. In heaving motion, there is no spatial distance variation, only an additional gravitational acceleration affecting the coolant flow. In rolling and pitching motions, the coolant flow in nuclear reactor systems is affected by an additional force in a non-inertial coordinate (Zhou et al., 2015). The additional force consists of three parts: tangential force, normal force and the Coriolis force, as shown in Fig. 10. Table 3 presents the comparison of heeling, heaving, rolling and pitching motions. It demonstrates that the heeling motion is the simplest, followed by the heaving motion. The rolling motion and pitching motion are the most complex.
For the viewpoint of thermal hydraulic considerations, it could be feasible to develop system analysis codes in ocean motions by modifying the existed codes. The modifying process includes two parts: theoretical models and field equations. Conventional models could be replaced with the models applicable in ocean motions. The use of CFD codes in complex channels and facilities could be explored to further explain natural circulation due to the small thermal driving head and flow instability in ocean motion due to thermal-induced instability and motion-induced instability. Understanding CHF in ocean motions is rather different when compared to the stationary state. Results obtained thus far indicated that the rolling CHF is lower than the stationary CHF. However, some experiments have shown that the rolling CHF increased in some thermal conditions. A more fundamental and extensive study is required to better understand CHF in ocean motions, particularly at high temperature and high pressure.

Thermal hydraulic considerations of the fourth generation nuclear reactor systems
The six different Generation IV reactor technologies that are currently being considered include: (i) Very High Temperature Reactor (VHTR), Super-Critical Water-cooled Reactor (SCWR), Lead-cooled Fast Reactor (LFR), Sodiumcooled Fast Reactor (SFR), Gas-cooled Fast Reactor (GFR) and Molten Salt Reactor (MSR). Goals of the Generation IV reactors are improved sustainability, economics, safety, and reliability together with proliferation resistance and physical protection. Sustainable energy generation involves meeting the clean air objectives, providing long term availability and effective utilization of fuel, and waste production should be minimized. In this paper, the thermal hydraulic challenges for SCWR, LFR and SFR are reviewed.
SCWR design is essentially an incremental evolution of the currently used Generation III/III and LWR. It can thus be treated as a combination of PWR and boiling water reactor (BWR). Within the SCWR pressure vessel, the working fluid does not change state as in a PWR, but rather the supercritical fluid/steam directly proceeds from rector pressure vessel to the turbine as for BWR. The utilisation of water above the critical point (T pc = 647 K, P cr = 22.064 MPa) can improve LWR efficiency by up to 30%. From the safety and operational perspective as well as thermal hydraulic considerations of the different SCWR designs from Japan, Canada, China, Russia, USA and Republic of Korea, the number of passes through the core influences the coolant, moderator and fuel behaviour. Although a one-pass core is a simpler and preferable design and provides natural circulation that will improve safety, there are three problems that are required to be resolved: (i) the axial peaking factor is large; (ii) it can be rather difficult to ensure a negative void reactivity coefficient at the end of the fuel life cycle; (iii) there is a huge increase in enthalpy along the channel that could lead to the overheating of fuel elements. Solutions to the above problems require the use of different fuel enrichment in the axial direction and/or additional moderation in the upper part of the core. A multi-pass flow in the core avoids negative void reactivity coefficients throughout the fuel life cycle as between each pass the flow mixes to avoid hot fuel rod spots from the high rise of enthalpy. The drawback of this concept is high pressure loss that reduces the effectiveness of natural circulation. Particularly, the safety analysis should ensure that the maximum cladding surface temperature is below 850 °C during abnormal transient conditions (Yamaji et al., 2006) and 1260 °C in the case of accidents that may lead to core damage (Okano et al., 1997) in order to avoid severe environmental repercussions such as radioactive material leakage. More details of the specific SCWR designs can be referred in Rowinski et al. (2018). Single-phase heat transfer plays a crucial role in the thermal-hydraulic evaluation of the nuclear reactor core of liquid metals fast reactor such as LFR. Design and safety analyses could be performed via system analysis and CFD codes. In these system analysis codes, the heat transport within the core can be accounted by means of correlations. However, the heat transfer correlations for liquid metals show a large spread as clearly demonstrated by Chandra et al. (2009). CFD codes can be used to evaluate possible local effects which cannot be derived from system analysis codes. Especially, the application of wire wraps commonly envisaged as spacer design for the fuel assemblies (Abderrahim et al., 2010) which requires a thorough knowledge of the local effects. The state-of-the-art for liquid metals fast reactor was shown to be LES simulations (see Fig. 6 for the different thermal hydraulics approaches within the CFD framework) for single pins up to a complete fuel assembly. LES simulations are currently used as reference data-set for the validation of computationally less demanding engineering resources such as RANS. The flow exiting the core is made up of the outlets of many different fuel assemblies. Liquid metal in these assemblies maybe heated up to different temperatures which could lead to temperature fluctuations on various above core structures. These temperature fluctuations may lead to thermal fatigue damage of the structures, an accurate characterization of the liquid metal flow structure in the above core region is thus very important. More details on the specific applications via system analysis and CFD codes on the fuel assembly and the pool thermal hydraulics being subdivided into the upper plenum and lower plenum can be found in Roelofs et al. (2013).
The prime candidate for large-scale implementation of breeder reactor technology is the SFR. For this type of Generation IV nuclear reactor system, liquid sodium is used as a coolant. Nevertheless, the heat transfer to liquid metals is significantly different from the heat transfer to water as liquid metals possess a very low Prandtl number. During the consideration of severe faults in this type of nuclear reactor system such as the unprotected loss of flow, loss of piping integrity, loss of heat sink, anticipated transient without scram and subassembly blockage, the boiling of sodium may appear, leading to dryout and even the melting of fuel bundles. Compared with the boiling of water in LWR or PWR technology, the boiling of liquid metal according to Sorokin et al. (1999) can be characterised in the following: the complex interaction of the internal factors in the system makes it difficult to accurately determine the incipient boiling superheat of liquid sodium under actual conditions; large vapor bubbles are formed in liquid sodium at several nucleation sites and the formation time of most vapor bubbles is within the waiting period; the growth of liquid sodium vapor bubbles can be explosive at a rate of about 10 m/s; the major two-phase flow patterns of liquid sodium are the same as that of conventional fluids and dispersed annular flow pattern dominates around the barometric pressure; the phase change of dispersed annular flow of liquid sodium in pipe is realized through the evaporation of liquid film rather than the formation of vapour bubbles (formed by boiling) on the wall surface; and the heat transfer coefficient can be very large.
For liquid sodium, the wall superheat of incipient boiling is higher than that of stable boiling by six times and about five times higher than that of water's incipient boiling due to: (1) the saturated vapor curve pressure gradient dp/dT of liquid sodium is relatively small within the universal pressure and temperature ranges; (2) the liquid sodium with active chemical property has strong self-cleaning effect on the heating surface, and furthermore the high invasion property of liquid sodium makes it difficult for the activation of surface nucleation; (3) the inert gases usually retain in wall-surface holes, thus favourable for vapor bubble nucleation, while hightemperature sodium tends to eliminate gas retention. Owing to its high heat conductivity and large vapourliquid specific volume difference (relative to water), the in-tube boiling process will mainly include annular flow and high vapour content, and critical boiling appears in the form of dryout. If the nuclear reactor systems suffer from flow loss, hot trap loss, transient overpower, transient underheating, local blockage or severe faults, the positive reaction effect of sodium bubbles due to boiling could make the fuel element melt or burn out. Specifically, the melting and burning damage of fuel element is closely related to the critical heat flux of liquid sodium (Subbotin et al., 1970). As a result of the special characteristics of sodium, different models for sodium boiling mechanism have been developed and meanwhile provided the corresponding solutions. Nevertheless, there remain many discrepancies and even conflicts between the proposed models. In summary, more studies are still required to be carried out for sodium twophase flow and heat transfer, especially in the aspect of sodium boiling heat transfer in bundle channels. An extensive review on the two-phase heat transfer and flow characteristics of liquid sodium can be found in Wu et al. (2018).

Summary
Nuclear safety analyses via the use of system analysis, subchannel and CFD codes to the present nuclear reactor technologies such as Generation II, III and III+ power reactors as well as future Generation IV reactor technologies are integral to the safe operation of nuclear reactor systems not only under normal conditions but also in the design of mitigating potential severe faults. Understanding the inherent feature of FBC that could lead to the occurrence of a range of flow instabilities in water-cooled and water-moderated nuclear reactor systems requires the critical understanding on the dynamics of different twophase flow regimes that exist in the primary and secondary cooling loops of Generation II, III and III+ reactor systems. A range of different experimental approaches have been performed to develop empirically determined flow regime maps through visualisation techniques and measure key parameters such as void fraction and interfacial area concentration through sophisticated methods in order to assess the validity of predictive two-phase flow models predominantly based on the interpenetrating media framework in system analysis, subchannel and CFD codes. Despite of the innumerable experiments and model developments that have been performed, a comprehensive description and the full prediction capability of CHF remains challenging and elusive. More work is still required to develop highly accurate empirical correlations and/or subchannel analysis as well as mechanistically based CFD models for design applications and for predicting transient behaviour of CHF. Potential use of nanofluids for CHF enhancement, small modular reactors on floating ocean platforms and Generation IV reactor technologies adopting coolants different from those of water-based nuclear reactor systems represents some of the key thermal hydraulic challenges being presented. During this exciting "trip to the moon" journey with regards to the many thermal hydraulic challenges, it can be concluded that there have certainly been important spinoffs that have been achieved despite the fact that the landing remains far reaching.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.