Modelling of Boiling Flows for Nuclear Thermal Hydraulics Applications—A Brief Review

: The boiling process is utterly fundamental to the design and safety of water-cooled ﬁssion reactors. Both boiling water reactors and pressurised water reactors use boiling under high-pressure subcooled liquid ﬂow conditions to achieve high surface heat ﬂuxes required for their operation. Liquid water is an excellent coolant, which is why water-cooled reactors can have such small sizes and high-power densities, yet also have relatively low component temperatures. Steam is in contrast a very poor coolant. A good understanding of how liquid water coolant turns into steam is correspondingly vital. This need is particularly pressing because heat transfer by water when it is only partially steam (‘nucleate boiling’ regime) is particularly e ﬀ ective, providing a great incentive to operate a plant in this regime. Computational modelling of boiling, using computational ﬂuid dynamics (CFD) simulation at the ‘component scale’ typical of nuclear subchannel analysis and at the scale of the single bubbles, is a core activity of current nuclear thermal hydraulics research. This paper gives an overview of recent literature on computational modelling of boiling. The knowledge and capabilities embodied in the surveyed literature entail theoretical, experimental and modelling work, and enabled the scientiﬁc community to improve its current understanding of the fundamental heat transfer phenomena in boiling ﬂuids and to develop more accurate tools for the prediction of two-phase cooling in nuclear systems. Data and insights gathered on the fundamental heat transfer processes associated with the behaviour of single bubbles enabled us to develop and apply more capable modelling tools for engineering simulation and to obtain reliable estimates of the heat transfer rates associated with the growth and departure of steam bubbles from heated surfaces. While results so far are promising, much work is still needed in terms of development of fundamental understanding of the physical processes and application of improved modelling capabilities to industrially relevant ﬂows.


Introduction
Understanding the boiling process is of great importance for the design and operation of light water reactors (LWRs).Both boiling water reactors (BWRs) and pressurised water reactors (PWRs) use nucleate boiling under high-pressure subcooled liquid flow conditions to achieve high surface heat fluxes required for their operation.A key requirement of these designs is to avoid the breakdown in heat transfer performance that can occur when a boiling surface becomes blanketed with vapour, potentially resulting in overheating, damage and release of fission products from damaged fuel rods.The breakdown in nucleate boiling heat transfer at high heat fluxes is referred to as critical heat flux (CHF), and it is imperative to be able to predict heat transfer in boiling flows to ensure normal reactor operation far from CHF conditions.A significant amount of work and effort are underway in the nuclear thermal hydraulics research community to develop simulation methods for flow boiling in LWR conditions and enable reliable analysis of boiling flows at the scale of nuclear reactor components.
The development of simulation techniques for boiling flows is at the core of current nuclear thermal hydraulics research.Current research trends entail an amalgamation of modelling-both computational and theoretical, and experimental work.The current review is focused on the modelling aspect.As we shall see, recent modelling research generated expertise and capabilities in at least two areas: (i) it enabled significant advancement in the development of practical engineering tools for the analysis of boiling flows and (ii) the parallel development and application of advanced physics-based modelling tools enabled gaining data and insights on fundamental processes in boiling fluids.
In this paper, recent literature on boiling phenomena and their modelling is reviewed and interpreted in the wider context of current nuclear thermal hydraulics research.Computational modelling of boiling is a buoyant area of research at the intersection of computational physics, thermal sciences and fluid mechanics and the state of the art is advancing rapidly.Although the topic has been reviewed as recently as 2014 [1], significant recent advancements prompted a new survey of the latest research developments.
The current paper is structured as follows: Section 2 provides an overview of the simulation methodology, embedded in the computational fluid dynamics (CFD) modelling framework, almost universally used for reactor analysis at the scale of the single plant components.In Section 3, the physical model forming the basis of all representations of boiling at a solid surface used in CFD simulation is discussed.In Section 4, a discussion of the accuracy of predictions of CFD simulation with basic wall-boiling modelling is presented.Attempts at developing more physically sound phenomenological models of boiling and their implementation in CFD simulation methods is discussed are reviewed.A more capable tool for the computational analysis of boiling phenomena, the so-called interface-capturing simulation approach, is discussed in Section 5, where results of the initial application of this recently developed methodology are reviewed.Outlook on future research needs to enable the development of nuclear systems relying on boiling is presented in Section 6.

Our Current Understanding of Vertical Upward Subcooled Flow Boiling
Computational modelling of boiling at the macroscopic scale typical of, for example, reactor sub-channels analyses is used to predict the hydrodynamic and thermal behaviour of non-equilibrium, diabatic liquid flows containing a population of steam bubbles.In these typically vertically oriented flows, the fluid enthalpy, initially below saturation ('subcooled' conditions) is increased by the addition of heat from the tube walls, which causes steam bubbles to be generated at the solid surface.A schematic of the flow is indicated in the diagram of Figure 1 [2], where z indicates the vertical distance along the channel, and where it is assumed that a vertical pipe is being heated by a constant and uniform heat flux.The physical situation sketched in Figure 1 allows identifying, in the lower section of a subchannel heated length, an 'onset of nucleate boiling' (ONB) point where conditions are reached at the surface for the formation of steam bubbles.In highly subcooled flows, these bubbles are short-lived and typically condense while still being attached to the solid surface.In these conditions, the boiling process does not cause appreciable production of voids.
Going up the heated length, the fluid enthalpy gradually increases due to the addition of heat at the pipe surface, until eventually a 'net vapour generation' (NVG) condition is reached whereby bubbles do not condense immediately after being generated, therefore contributing to the creation of permanent voids in the fluid.The highly subcooled flow typical of the section between the ONB and NVG point is characterised by a negligible presence of voids, typically concentrated at the heated surface, and by almost extreme processes of evaporation from the superheated liquid to the vapour near the wall, and of condensation from the vapour to the subcooled liquid at a short distance from the wall.During this flow regime, bubbles condense at a short distance from the solid surface, or sometimes even without lifting off the surface, whereas above the NVG point, bubbles penetrate some distance into the bulk liquid and therefore contribute to the generation of a permanent vapour phase.some distance into the bulk liquid and therefore contribute to the generation of a permanent vapour phase.
Eventually, above a certain height, the entirety of the liquid phase is at or above the saturation temperature corresponding to the prevailing pressure.The corresponding saturated boiling regime, characterised by a significant presence of voids, can be thought of as representative of the upper section of a PWR hot channel or the lower-mid section of a BWR channel.Conditions of interest for nuclear applications are characterised by high pressures (up to 160 bar), large liquid flow velocities of a few meters per second, and typical bubble departure sizes of order 0.01-0.1 mm.Boiling occurs in pipes of about 10 mm of diameters and up to a few metres of length, and bubbles are generated at the inner surface of the pipe at an areal density of hundreds or perhaps thousands of bubbles per square centimetre of pipe wall, at a frequency of hundreds of bubbles per second.The flow is also highly turbulent, and the temperature variation in any pipe cross-section is of several tens of degrees.The daunting complexity of this physical picture inevitably results in a series of simplifications in its physical modelling, which are required to formulate tractable and numerically solvable computational models of the process.

Overview of a Practical Simulation Method for Component-Scale Analysis of Boiling
The main approximation involved in the formulation of practical component scale models of boiling flows is to assume that the liquid and vapour phases can be represented by continuous interpenetrating fields.This is an inevitable consequence of the mismatch between typical bubble sizes (smaller than a fraction of a millimetre) and the dimensions of the thermal hydraulic systems of interest (pipes several metres long and a few centimetres in diameter).In such conditions, tracking the behaviour of the single bubbles is not practically feasible, and the bubbles are replaced by a representation of the void fraction associated with their presence in the fluid.The approach is usually termed a 'Eulerian-Eulerian' [3,4], or two-fluid, method.
Two-fluid methods model spatially and temporally averaged quantities to represent thermal and hydrodynamic characteristics of two-phase systems.The fluid mixture behaviour is modelled using one set of conservation equations (for mass, momentum and energy) for each phase [5].For boiling flows, one set is used to model the vapor (i.e., the bubbles) phase, while the other set is solved The bottom panel shows a diagram of the likely distribution of steam bubbles generated at a heated wall in vertical subcooled boiling.In the top panel, the likely distribution of cross-section averaged void fraction is shown.Flow regimes are identified depending on void fraction values: 'single-phase convection'-zero void fraction; 'highly subcooled'-bubbles are present only near the wall and condense as soon as released into the flow; 'low subcooled'-permanence of steam bubbles in a body of partly subcooled liquid; 'saturated boiling'-the fluid is entirely at or above the saturation temperature and bubbles fill up the pipe cross-section.
Eventually, above a certain height, the entirety of the liquid phase is at or above the saturation temperature corresponding to the prevailing pressure.The corresponding saturated boiling regime, characterised by a significant presence of voids, can be thought of as representative of the upper section of a PWR hot channel or the lower-mid section of a BWR channel.
Conditions of interest for nuclear applications are characterised by high pressures (up to 160 bar), large liquid flow velocities of a few meters per second, and typical bubble departure sizes of order 0.01-0.1 mm.Boiling occurs in pipes of about 10 mm of diameters and up to a few metres of length, and bubbles are generated at the inner surface of the pipe at an areal density of hundreds or perhaps thousands of bubbles per square centimetre of pipe wall, at a frequency of hundreds of bubbles per second.The flow is also highly turbulent, and the temperature variation in any pipe cross-section is of several tens of degrees.The daunting complexity of this physical picture inevitably results in a series of simplifications in its physical modelling, which are required to formulate tractable and numerically solvable computational models of the process.

Overview of a Practical Simulation Method for Component-Scale Analysis of Boiling
The main approximation involved in the formulation of practical component scale models of boiling flows is to assume that the liquid and vapour phases can be represented by continuous interpenetrating fields.This is an inevitable consequence of the mismatch between typical bubble sizes (smaller than a fraction of a millimetre) and the dimensions of the thermal hydraulic systems of interest (pipes several metres long and a few centimetres in diameter).In such conditions, tracking the behaviour of the single bubbles is not practically feasible, and the bubbles are replaced by a representation of the void fraction associated with their presence in the fluid.The approach is usually termed a 'Eulerian-Eulerian' [3,4], or two-fluid, method.
Two-fluid methods model spatially and temporally averaged quantities to represent thermal and hydrodynamic characteristics of two-phase systems.The fluid mixture behaviour is modelled using one set of conservation equations (for mass, momentum and energy) for each phase [5].For boiling flows, one set is used to model the vapor (i.e., the bubbles) phase, while the other set is solved for the liquid phase.The unique aspect of the approach is that it enables practical, albeit approximate, modelling of interaction between the phases.These interfacial models represent the mass, momentum and energy transfers through the interface between the phases.The main limitation of the approach is that it requires specification of closure laws for interfacial mass, momentum and energy transfer, and for interfacial area transport.These interfacial interactions are of crucial importance as they determine the rate of phase change, and the degree of mechanical and thermal non-equilibrium between phases [5].
Mass transfer between phases is present in the form of evaporation due to the addition of heat from a solid surface and due to vapor bubbles being condensed in the bulk liquid (which typically is at a temperature below saturation).From a hydrodynamic point of view, important interfacial effects between vapour (disperse phase, in the form of bubbles) and liquid (continuous phase) are due to vapour-liquid drag force [6] as well as other possible interactions such as lift [7] and wall lubrication [8] forces, and turbulence-assisted bubble dispersion [9].Furthermore, the distribution of local bubble sizes in the subcooled liquid flow is expected to be strongly influenced by the complex bubble behaviours in the two-phase flows occurring near the heated wall [10].
A necessary simplification of the two-fluid approach is dictated by the highly turbulent nature of the flows.The time step and mesh requirements for accurate description of the turbulent flow behaviour with scale-resolving simulation (using for example direct numerical or large eddy simulation) would not allow repetition in reasonable time of simulation runs for different design parameters, a practice required by engineering analysis.Therefore, the two-fluid approach is normally coupled to a turbulence model capturing the temporally averaged behaviour of the flow, using the Reynolds-averaged formulation of the Navier-Stokes equations ('RANS' approach) [11].
For boiling flows of practical interest, additional special care is required by the modelling of wall-to-fluid energy transfer in the presence of bubbles [12], which is described in the following section.

Current Understanding of the Cycle of Processes Associated with Boiling at a Surface
A general first-principles based mathematical model of boiling does not exist yet due to the complex multi-physics and multi-scale nature of the problem [13].Our current physical understanding of boiling at a surface is limited and mainly derived from laboratory experiment in conditions quite distant from real applications.With this in mind, it is however possible to sketch a cycle of processes associated with the generation and departure of steam bubbles at a solid surface, as shown in Figure 2a.It is believed that bubbles are created at locations ('nucleation sites') on the heated surface where crevices or impurities create the conditions for the permanence of microscopic voids.Bubbles grow at these nucleation sites due to evaporation at the bubble curved surface ('superheated layer' evaporation).In low-pressure boiling, a thin liquid film ('microlayer) is often present beneath the bubble [14], which evaporates due to the addition of heat directly from the wall, thus contributing as well to bubble formation.Evaporation causes the bubble to grow in size until eventually buoyancy forces are sufficient to lift the bubble off the surface.In the process, cold fluid is brought in the vicinity of the wall, causing the local temperature to fall below the value typically required to sustain bubble growth.During this quiescent period, the addition of heat from the wall causes the local wall and liquid temperature to increase around the nucleation site, until conditions are again met for the growth of another bubble.

Basic Wall-Boiling Model for Eulerian-Eulerian Simulation
Engineering applications require us to compute the heat transfer rates due to cooling via a boiling fluid, specifically for vertical bubbly flow at low quality with generation of bubbles near a solid surface ('subcooled flow boiling'), with little or no net generation of vapour [15,16].

Basic Wall-Boiling Model for Eulerian-Eulerian Simulation
Engineering applications require us to compute the heat transfer rates due to cooling via a boiling fluid, specifically for vertical bubbly flow at low quality with generation of bubbles near a solid surface ('subcooled flow boiling'), with little or no net generation of vapour [15,16].Computational codes such as STAR-CCM+ [4] and ANSYS Fluent [17] for engineering analysis of boiling flows at component scale, typically used in nuclear fuel channel analysis, are augmented with subgrid scale models accounting for the heat transfer modes associated with the formation of steam bubbles at the heated wall.The approach that forms the basis of these treatments is the so-called 'heat flux partitioning' method of Kurul and Podowski [18], also known as the 'Rensselaer Polytechnic Institute (RPI)' model.The physical basis of the model, sketched in Figure 2b, postulates two additional heat transfer modes typical of boiling conditions, which augment the ordinary single-phase convective heat transfer to the liquid.These additional heat transfer modes are associated with the formation of bubbles, an 'evaporative heat flux' representing the latent heat transfer due to bubble generation, and to the process of bubble departure, a 'quench' heat flux associated with the departure of a bubble causing cold fluid to make contact with the heated wall.The basic assumption of the model is to treat these heat transfer modes separately, linking the magnitude of the three heat transfer components to the areal density of bubble nucleation sites, departure diameter, and lift-off frequency of steam bubbles.The model requires input data on these parameters, typically taken from an experiment [19].

Manual Assessment of the RPI Model
Application of the RPI model requires knowledge of various fundamental aspects of the boiling process, typically provided in the form of bubble departure diameters, nucleation site densities and bubble departure frequencies.First-principles based modelling tools to compute these quantities are still in the early stages of development (as discussed in Section 5), therefore, it is customary to use experimental correlations in order to provide input data required by RPI.
When compared against detailed experimental observations of the bubble formation process, experimental correlations, widely used in RPI-based component scale analyses, usually return poor predictions.Figure 3, comparing bubble departure diameters measured by Goel et al. [20] against predictions of the Kocamustafaogullari [21] correlation, provides a typical example of the magnitude of the discrepancy [20].
However, a manual evaluation [15] of RPI-modelled heat transfer showed that errors due to inaccuracy of correlations, employed to supply bubble departure diameters, frequencies and nucleation site densities, are likely cancel out, and that although accuracy of the raw input data from correlations is usually low, derived quantities such as the wall temperature can adequately be predicted by the RPI model.This outcome is indicated in Figure 4, adapted from Thakrar et al. [15], showing a comparison between measured wall temperatures, their values predicted using manual evaluation of the RPI model and existing boiling heat transfer correlations for a typical boiling validation case.The blue crosses in the figure are suggestive of good agreement between RPI-modelled and measured wall temperatures, confirming the physical basis of the RPI approach.

Importance of the Wall-Boiling Model
Application of RPI requires knowledge of various fundamental aspects of the boiling process to work, typically provided via the bubble departure diameter, nucleation site density and frequency.Modelling tools to compute these quantities are still in the early stages of development, therefore, it is customary to use correlations in order to provide input data required by RPI.Typical CFD

Importance of the Wall-Boiling Model
Application of RPI requires knowledge of various fundamental aspects of the boiling process to work, typically provided via the bubble departure diameter, nucleation site density and frequency.Modelling tools to compute these quantities are still in the early stages of development, therefore, it is customary to use correlations in order to provide input data required by RPI.Typical CFD

Importance of the Wall-Boiling Model
Application of RPI requires knowledge of various fundamental aspects of the boiling process to work, typically provided via the bubble departure diameter, nucleation site density and frequency.Modelling tools to compute these quantities are still in the early stages of development, therefore, it is customary to use correlations in order to provide input data required by RPI.Typical CFD simulation results, as shown in Figure 5 from Colombo et al. [4], indicate that predictions of important quantities such as the volume fraction distribution can vary to a large extent depending on the particular treatment (in Figure 5, different bubble departure correlations) used for generating input data to the RPI wall-boiling model.simulation results, as shown in Figure 5 from Colombo et al. [4], indicate that predictions of important quantities such as the volume fraction distribution can vary to a large extent depending on the particular treatment (in Figure 5, different bubble departure correlations) used for generating input data to the RPI wall-boiling model.If the various adjustable coefficients, most importantly bubble departure diameter values, are carefully selected, the approach enables simulation in challenging geometries and industrially relevant conditions.For example, boiling in unusual geometries, such as non-circular channels [24], which are still challenging from a hydrodynamic point of view, can be reliably predicted, within the inevitable limitations of a largely semi-empirical approach, which requires careful tuning of input parameters in order to obtain reliable estimates of the various quantities of interest, such as void fraction distributions and boiling heat transfer rates.
In order to improve accuracy and generality of CFD predictions using the RPI model, it is of crucial importance to provide the RPI wall-boiling model with accurate input data.To this end, efforts are focused on improving accuracy of input data for the bubble departure diameter, which is the single most sensitive parameter of the RPI model.

Development of Wall-Boiling Models
CFD simulation of boiling using RPI is sensitive to input data provided to RPI heat transfer submodels.These depend strongly on input bubble departure diameter value, for which experimental databases are scarce and do not usually cover the range of parameter space typical of industrially relevant flow conditions.If the various adjustable coefficients, most importantly bubble departure diameter values, are carefully selected, the approach enables simulation in challenging geometries and industrially relevant conditions.For example, boiling in unusual geometries, such as non-circular channels [24], which are still challenging from a hydrodynamic point of view, can be reliably predicted, within the inevitable limitations of a largely semi-empirical approach, which requires careful tuning of input parameters in order to obtain reliable estimates of the various quantities of interest, such as void fraction distributions and boiling heat transfer rates.
In order to improve accuracy and generality of CFD predictions using the RPI model, it is of crucial importance to provide the RPI wall-boiling model with accurate input data.To this end, efforts are focused on improving accuracy of input data for the bubble departure diameter, which is the single most sensitive parameter of the RPI model.

Development of Wall-Boiling Models
CFD simulation of boiling using RPI is sensitive to input data provided to RPI heat transfer submodels.These depend strongly on input bubble departure diameter value, for which experimental No mechanistic model of the bubble nucleation and growth process exists, and, therefore, it is not possible to predict bubble of formation and release in any given conditions; however, a number of semi-mechanistic phenomenological approaches have been developed to compute the diameters of steam bubbles at departure from a heated surface in flow boiling.
One of these approaches is the popular 'force-balance' method of Klausner [25][26][27][28][29], which, as indicated in Figure 6, computes the size of a bubble at departure from a solid surface based on evaluation of hydrodynamic (drag force F DF , lift force F L and liquid reaction to bubble expansion F H ), surface tension, gravity and wall-adhesion forces.
A different approach, based on thermodynamic considerations, [30][31][32][33], allows tracking the dynamic contact angle at the base of a bubble, and therefore to compute the temporally varying radial extent of the bubble base diameter and of the surface tension force, which can be crucially important for very small bubbles (likely as small as perhaps 10 micron in typical PWR conditions).Figure 7 shows an assessment of various bubble departure models, showing generally large errors and therefore a perceived need to elucidate further the process of bubble formation and departure.
Recent efforts [27,35] aimed to improve CFD prediction of boiling heat transfer using data generated with semi-mechanistic bubble departure models.In [24] in particular, popular 'force balance' models were used to compute bubble departure diameters and used in CFD simulation of boiling flows.The models were assessed via comparison with established validation cases close to real PWR conditions, showing that approach can return predictions in good agreement with experimental trends Figure 8. Massive variation can be observed in values of departure diameters as predicted by the various phenomenological models or correlations and comparatively much smaller differences in 'derived' quantities such as the wall temperature computed by the CFD model, although significant error can arise depending on the choice of input values for the bubble departure diameter.No mechanistic model of the bubble nucleation and growth process exists, and, therefore, it is not possible to predict bubble of formation and release in any given conditions; however, a number of semi-mechanistic phenomenological approaches have been developed to compute the diameters of steam bubbles at departure from a heated surface in flow boiling.
One of these approaches is the popular 'force-balance' method of Klausner [25][26][27][28][29], which, as indicated in Figure 6, computes the size of a bubble at departure from a solid surface based on evaluation of hydrodynamic (drag force , lift force and liquid reaction to bubble expansion ), surface tension, gravity and wall-adhesion forces.
A different approach, based on thermodynamic considerations, [30][31][32][33], allows tracking the dynamic contact angle at the base of a bubble, and therefore to compute the temporally varying radial extent of the bubble base diameter and of the surface tension force, which can be crucially important for very small bubbles (likely as small as perhaps 10 micron in typical PWR conditions).Figure 7 shows an assessment of various bubble departure models, showing generally large errors and therefore a perceived need to elucidate further the process of bubble formation and departure.[32] of the energy based model of Ardron et al. [30].The black lines indicate results of application of the Klausner [25] and Yun [34] models, and of the empirical fit to the Klausner model by Sugrue et al. [28].From [32].
Recent efforts [27,35] aimed to improve CFD prediction of boiling heat transfer using data generated with semi-mechanistic bubble departure models.In [24] in particular, popular 'force balance' models were used to compute bubble departure diameters and used in CFD simulation of boiling flows.The models were assessed via comparison with established validation cases close to real PWR conditions, showing that approach can return predictions in good agreement with experimental trends Figure 8. Massive variation can be observed in values of departure diameters as predicted by the various phenomenological models or correlations and comparatively much smaller differences in 'derived' quantities such as the wall temperature computed by the CFD model, although significant error can arise depending on the choice of input values for the bubble departure diameter.

The Interface-Capturing Simulation Approach
Interface capturing methods belong to a class of computational techniques devised to track the evolution of fluid interfaces, typically coupled to the CFD representation of fluid flow in order to predict two-phase fluid processes [36].This class of methods is a particularly promising candidate for enabling ab initio prediction of two-phase flows owing to a number of desirable properties.With the interface capturing approach, the behaviour of the vapour-liquid interface is modelled via solution of the fundamental two-phase flow equations in their instantaneous formulation, which enables modelling arbitrary interface configurations and flow regimes.The effect of surface tension is modelled locally at the interface, which provides a framework for application of the methodology to small-scale flows in which surface tension forces play a dominant role.Finally, and perhaps most importantly, the interface capturing method enables straightforward extension to modelling heat transfer in both fluid phases and the computation of interphase mass transfer locally at the interface, which is crucial for boiling phenomena.
Methods vary depending on the particular technique employed to mark the different fluid

The Interface-Capturing Simulation Approach
Interface capturing methods belong to a class of computational techniques devised to track the evolution of fluid interfaces, typically coupled to the CFD representation of fluid flow in order to predict two-phase fluid processes [36].This class of methods is a particularly promising candidate for enabling ab initio prediction of two-phase flows owing to a number of desirable properties.With the interface capturing approach, the behaviour of the vapour-liquid interface is modelled via solution of the fundamental two-phase flow equations in their instantaneous formulation, which enables modelling arbitrary interface configurations and flow regimes.The effect of surface tension is modelled locally at the interface, which provides a framework for application of the methodology to small-scale flows in which surface tension forces play a dominant role.Finally, and perhaps most importantly, the interface capturing method enables straightforward extension to modelling heat transfer in both fluid phases and the computation of interphase mass transfer locally at the interface, which is crucial for boiling phenomena.
Methods vary depending on the particular technique employed to mark the different fluid phases on the computational grid used by a CFD modelling technique.Differences have been addressed extensively elsewhere (see, for example, books by Tryggvason et al. [36] and the review of Nandi et al. [3]), and, for present purposes three established methods are distinguished: -With the volume of fluid (VOF) method [37], the volume fraction of the 'primary' fluid is used to distinguish the two phases; - The level set (LS) method [38] identifies the interface as the zero level of a function representing the shortest distance from the interface; - The front tracking (FT) [39] method describes the interface as a set of massless particles moved around by the fluid velocity field.
The power of the interface capturing approach lies in its ability to model microscopic interfacial processes locally, which enables predicting from first principles the overall flow behaviour, typically involving length scales of several millimetres.This aspect is recognized to be particularly important for the modelling of boiling [40] phenomena.In particular, interface capturing methods provide a framework for modelling thermally driven interfacial mass transfer phenomena [41,42], a typical case being the growth of steam bubbles due to the addition of heat from a solid surface.Understanding these processes is fundamental to the modelling of boiling phenomena, therefore, the interface capturing approach is very promising for application to modelling of boiling in industrial configurations.
To date reliable predictions have been possible with interface capturing methods only for ideal test cases involving one or at most a few bubbles growing on a heated surface, as indicated in Figure 9. Application of these methods is computationally intensive, requiring massive computational resources, which so far has been the single most important factor limiting applicability of the interface capturing method.As a consequence, the method has to date been applied only in laboratory conditions, and not yet in industrially relevant flow configurations.-With the volume of fluid (VOF) method [37], the volume fraction of the 'primary' fluid is used to distinguish the two phases; - The level set (LS) method [38] identifies the interface as the zero level of a function representing the shortest distance from the interface; - The front tracking (FT) [39] method describes the interface as a set of massless particles moved around by the fluid velocity field.
The power of the interface capturing approach lies in its ability to model microscopic interfacial processes locally, which enables predicting from first principles the overall flow behaviour, typically involving length scales of several millimetres.This aspect is recognized to be particularly important for the modelling of boiling [40] phenomena.In particular, interface capturing methods provide a framework for modelling thermally driven interfacial mass transfer phenomena [41,42], a typical case being the growth of steam bubbles due to the addition of heat from a solid surface.Understanding these processes is fundamental to the modelling of boiling phenomena, therefore, the interface capturing approach is very promising for application to modelling of boiling in industrial configurations.
To date reliable predictions have been possible with interface capturing methods only for ideal test cases involving one or at most a few bubbles growing on a heated surface, as indicated in Figure 9. Application of these methods is computationally intensive, requiring massive computational resources, which so far has been the single most important factor limiting applicability of the interface capturing method.As a consequence, the method has to date been applied only in laboratory conditions, and not yet in industrially relevant flow configurations.

Further Developments for Extension to Boiling Conditions
Interface capturing methods were originally devised for isothermal flows with no mass transfer and in the absence of solid boundaries.Therefore, further developments in the methodology are required to model strong interphase mass transfer and the effect of solid boundaries, which are typical aspects of boiling flows.

Modelling Mass Transfer
Various approaches are possible [44] for capturing mass transfer at a fluid interface.For the case of evaporation mass transfer typical of boiling flows, essentially two classes of approaches have been developed, used in conjunction with various interface capturing techniques, as summarised in Table 1.One class of methods assumes thermally driven mass transfer due to the conduction of heat towards (for evaporation) or away from (for condensation) the vapour-liquid interface.These

Further Developments for Extension to Boiling Conditions
Interface capturing methods were originally devised for isothermal flows with no mass transfer and in the absence of solid boundaries.Therefore, further developments in the methodology are required to model strong interphase mass transfer and the effect of solid boundaries, which are typical aspects of boiling flows.

Modelling Mass Transfer
Various approaches are possible [44] for capturing mass transfer at a fluid interface.For the case of evaporation mass transfer typical of boiling flows, essentially two classes of approaches have been developed, used in conjunction with various interface capturing techniques, as summarised in Table 1.One class of methods assumes thermally driven mass transfer due to the conduction of heat towards (for evaporation) or away from (for condensation) the vapour-liquid interface.These methods are listed in Table 1 as 'conduction in interface cells' methods, whereby the conductive heat flux driving mass transfer is approximated as continuous variation across the computational cells containing the interface, and 'heat flux balance' methods, relying on an explicit reconstruction of the generally discontinuous heat fluxes on both sides of the interface.With a different approach, the rate of mass transfer is evaluated using an approximation of the net interphase mass flux derived from basic kinetic theory ('kinetic model').Various implementations of the method have been proposed [45][46][47] based on the original interface capturing method of Hardt and Wondra [48].The 'kinetic' and 'conduction' models are the most popular approaches for realistic simulation of boiling and have been used in various combinations with different approaches for capturing interface behaviour.Not as well established, although very powerful as demonstrated via comparison with established benchmark cases [49], is a radically different method for capturing mass transfer, indicated as the 'asymptotic relaxation model' in Table 1, and derived from the phase field representation of two-phase flows.Merits and shortcomings of the various approaches for modelling mass transfer depend on flow conditions, and a single framework that is general enough to be applied to any boiling flow has not been developed yet.

Modelling Bubble-Wall Interaction
In boiling flows of practical interest, steam bubbles are generated at a solid surface due to the addition of heat from the wall.Therefore wettability [61,62], or 'contact angle' [30], of the surface plays an important role in bubble behaviour.On the one hand, contact angles were found to determine the size of steam bubbles at departure from the surface [30,32].On the other hand, contact angle values determine the likelihood of the formation of liquid 'microlayers' [63,64] beneath growing steam bubbles.Evaporation of these microlayers contributes significantly to the generation of vapour in a bubble, and their hydrodynamic and thermal behaviour is still being researched.
Another important area of active research is focussed on the evaporation process as it happens near the solid surface.Two evaporation mechanisms have been so far identified, one taking place at the solid-vapour-liquid contact line [65] and involving the effect of molecular interactions between the solid substrate and the evaporating vapour-liquid interface, and one in which a liquid microlayer [66][67][68] of a few microns in thickness is observed to form and subsequently deplete due to evaporation on the heated surface.These are still poorly understood phenomena currently being investigated from a fundamental perspective.

Application of Interface-Capturing Simulation
Application of microscopic bubble simulation with interface capturing methods has to date been limited to modelling aspects of the pool boiling phenomenon in laboratory conditions.Extension to the more challenging flow boiling conditions is at present outside the remit of current modelling techniques.Nevertheless, recently reported applications [67,69] enabled gaining unprecedented insight on the hydrodynamics and heat transfer of bubble behaviour in boiling.

Computation of Bubble Departure Diameters and Frequencies
Interface-capturing simulation has been successfully used for mechanistic prediction of bubble growth [30] and departure [68,70], with the data thus generated available for extracting bubble departure diameters and frequencies directly from simulation [71], enabling their subsequent use as input data for RPI-based component scale models.Mechanistic models of the single bubbles have been devised to capture the growth of a bubble from an assumed 'seed' of vapour located at the heat transfer surface.The bubble nucleation and early formation processes are outside the remit of these methods, which need to be augmented via coupling with molecular mechanics [72] simulation in order to develop self-consistent models of the nucleation cycle.Nevertheless, use of nucleate boiling simulation from a pre-set seed of vapour enabled prediction and experimental validation of bubble growth in pressurized water close to typical reactor operation conditions [71].Comparison with experiment and simulation at atmospheric pressure [68] confirmed that the bubble departure diameter decreases dramatically with increasing system pressure, a fact with far-reaching implications on the expected thermal behaviour of reactor systems, which rely on high-pressure boiling of water.To get a sense of scales, according to our current understanding of bubble departure diameters as derived from interface-capturing simulations, bubble diameters of approximately 10 microns are expected in PWR conditions, whereas in atmospheric conditions, typical bubble sizes are of a few millimetres.

Surface Phenomena
Interface-capturing simulation augmented with conjugate heat transfer to model the thermal response of a solid substrate has been applied to analyse aspects of the boiling process influenced by surface variables.A typical example is application to modelling the process of formation and depletion of a liquid microlayer beneath a growing bubble, as shown in Figure 10.Hänsch and Walker [67,73] developed a hydrodynamic, level set-based model of bubble behaviour and coupled it to a simplified thermal model of the microlayer beneath the bubble.The coupled simulation was used to track the depletion of the film and to model transient heat conduction in the solid and its thermal response to the latent heat transfer associated with film depletion, which indicated that an as yet unknown thermal resistance phenomenon associated with the evaporation process [66] is the main factor limiting the rate of microlayer depletion.Another application of conjugate heat transfer modelling, shown in Figure 11, is that of Giustini et al. [69,74], who applied mechanistic single bubble simulation to study the thermal response of the solid surface to the bringing of cold fluid in the vicinity of the wall caused by the bubble departure process ('quench heat transfer' [11]).Comparison with direct measurements of the spatial and temporal distribution of the heat flows beneath the bubble confirmed the magnitude of quench heat transfer rates due to the thermal response of a hot solid surface in contact with cold liquid in the wake of a departing bubble.Another application of conjugate heat transfer modelling, shown in Figure 11, is that of Giustini et al. [69,74], who applied mechanistic single bubble simulation to study the thermal response of the solid surface to the bringing of cold fluid in the vicinity of the wall caused by the bubble departure process ('quench heat transfer' [11]).Comparison with direct measurements of the spatial and temporal distribution of the heat flows beneath the bubble confirmed the magnitude of quench heat transfer rates due to the thermal response of a hot solid surface in contact with cold liquid in the wake of a departing bubble.The key requirement of all nuclear reactor designs is to avoid the breakdown in heat transfer performance that occurs when the surface of nuclear fuel rods becomes blanketed by a film of vapour.As noted, the phenomenon known as critical heat flux (CHF) is the single most important design constraint to the development of nuclear systems and of light water reactors (LWRs) in particular.The ultimate aim of boiling simulation is prediction of the occurrence of CHF phenomena and their The key requirement of all nuclear reactor designs is to avoid the breakdown in heat transfer performance that occurs when the surface of nuclear fuel rods becomes blanketed by a film of vapour.As noted, the phenomenon known as critical heat flux (CHF) is the single most important design constraint to the development of nuclear systems and of light water reactors (LWRs) in particular.The ultimate aim of boiling simulation is prediction of the occurrence of CHF phenomena and their dependence on thermal hydraulics and surface-related parameters.The behaviour of boiling phenomena as the CHF limit is approached is still poorly understood, a fact that is indeed the main drive for current research efforts on the modelling of boiling.In this context, interface-capturing simulation has been demonstrated to be a very promising tool for mechanistic prediction of CHF.To this end, Sato and Niceno [56,75] developed a simulation tool for modelling the collective behaviour of a population of steam bubbles seeded at pre-set nucleation sites on a heated substrate, and used it to predict the occurrence of CHF via vapour blanketing of the surface in pool boiling laboratory conditions.Shown in Figure 12, this is the only example reported to date of mechanistic modelling of the collective behaviour of steam bubbles as the CHF limit is approached, and the method holds great promise for application to more challenging flow configurations typical of reactor operations.
this end, Sato and Niceno [56,75] developed a simulation tool for modelling the collective behaviour of a population of steam bubbles seeded at pre-set nucleation sites on a heated substrate, and used it to predict the occurrence of CHF via vapour blanketing of the surface in pool boiling laboratory conditions.Shown in Figure 12, this is the only example reported to date of mechanistic modelling of the collective behaviour of steam bubbles as the CHF limit is approached, and the method holds great promise for application to more challenging flow configurations typical of reactor operations./ , tracking the evolution of a few bubbles for two seconds until the bubbles merge and form a vapour blanket on the solid surface.

Outlook-Future Issues
1) The two-fluid approach to modelling flow boiling presented in this review relies on a number of approximations and empirical parameters that limit the applicability of the approach, which should be considered obsolete, as more capable and physically consistent methods have now been developed.2) Application of modern interface capturing methods to problems typical of boiling in laboratory conditions, typically in low-pressure, low-subcooling pool boiling mode, enabled gaining unprecedented insight on the process of bubble formation and release at a surface.However, none of the interface capturing methods discussed in this review are yet applicable to highly turbulent subcooled bubbly flows typical of reactor operations.3) Extensions to modelling the behaviour of the solid-liquid-vapour contact line at the base of a steam bubble are in their early stages of development, and a general model that is applicable to any fluid or surface material does not yet exist.4) Modelling of evaporation at the phase boundary (e.g., the curved surface of a bubble) has to date been possible only for the case of thermally driven evaporation in near-equilibrium conditions.Thus, efforts should be pursued to extend current modelling capabilities to nonequilibrium conditions and cases where both dynamic and thermal effects determine bubble behaviour.
Funding: G.G. is supported by the United Kingdom Engineering and Physical Sciences Research Council (EPSRC) through grant EP/T027061/1.

Conflicts of Interest:
The author declares no conflict of interest.
Figure 12.Simulation of the collective behaviour of steam bubbles in pool boiling [56] at a heat flux of 300 KW/m 2 , tracking the evolution of a few bubbles for two seconds until the bubbles merge and form a vapour blanket on the solid surface.

Outlook-Future Issues
(1) The two-fluid approach to modelling flow boiling presented in this review relies on a number of approximations and empirical parameters that limit the applicability of the approach, which should be considered obsolete, as more capable and physically consistent methods have now been developed.(2) Application of modern interface capturing methods to problems typical of boiling in laboratory conditions, typically in low-pressure, low-subcooling pool boiling mode, enabled gaining unprecedented insight on the process of bubble formation and release at a surface.However, none of the interface capturing methods discussed in this review are yet applicable to highly turbulent subcooled bubbly flows typical of reactor operations.(3) Extensions to modelling the behaviour of the solid-liquid-vapour contact line at the base of a steam bubble are in their early stages of development, and a general model that is applicable to any fluid or surface material does not yet exist.(4) Modelling of evaporation at the phase boundary (e.g., the curved surface of a bubble) has to date been possible only for the case of thermally driven evaporation in near-equilibrium conditions.Thus, efforts should be pursued to extend current modelling capabilities to non-equilibrium conditions and cases where both dynamic and thermal effects determine bubble behaviour.

Figure 1 .
Figure 1.Our current understanding of nucleate flow boiling in nuclear reactor conditions[2].The bottom panel shows a diagram of the likely distribution of steam bubbles generated at a heated wall in vertical subcooled boiling.In the top panel, the likely distribution of cross-section averaged void fraction is shown.Flow regimes are identified depending on void fraction values: 'single-phase convection'-zero void fraction; 'highly subcooled'-bubbles are present only near the wall and condense as soon as released into the flow; 'low subcooled'-permanence of steam bubbles in a body of partly subcooled liquid; 'saturated boiling'-the fluid is entirely at or above the saturation temperature and bubbles fill up the pipe cross-section.

Figure 1 .
Figure 1.Our current understanding of nucleate flow boiling in nuclear reactor conditions [2].The bottom panel shows a diagram of the likely distribution of steam bubbles generated at a heated wall in vertical subcooled boiling.In the top panel, the likely distribution of cross-section averaged void fraction is shown.Flow regimes are identified depending on void fraction values: 'single-phase convection'-zero void fraction; 'highly subcooled'-bubbles are present only near the wall and condense as soon as released into the flow; 'low subcooled'-permanence of steam bubbles in a body of partly subcooled liquid; 'saturated boiling'-the fluid is entirely at or above the saturation temperature and bubbles fill up the pipe cross-section.

Inventions 2020, 5 , 19 Figure 2 .
Figure 2. (a) Current understanding of the cycle of processes associated with boiling; (b) basic representation of the main mechanisms of wall heat transfer.From Thakrar et al. [15].

Figure 2 .
Figure 2. (a) Current understanding of the cycle of processes associated with boiling; (b) basic representation of the main mechanisms of wall heat transfer.From Thakrar et al. [15].

Figure 3 .
Figure 3.Comparison between bubble departure diameters generated with the popular correlation of Kocamustafaogullari[21], widely used in Rensselaer Polytechnic Institute (RPI) modelling of wall boiling, and measured bubble departure diameters in conditions of the work in[20].

Figure 4 .
Figure 4. Comparison between measured wall temperatures, their values predicted using manual evaluation of the RPI model and boiling heat transfer correlations, for a typical boiling validation case.From[15].

Figure 3 .
Figure 3.Comparison between bubble departure diameters generated with the popular correlation of Kocamustafaogullari[21], widely used in Rensselaer Polytechnic Institute (RPI) modelling of wall boiling, and measured bubble departure diameters in conditions of the work in[20].

Figure 3 .
Figure 3.Comparison between bubble departure diameters generated with the popular correlation of Kocamustafaogullari[21], widely used in Rensselaer Polytechnic Institute (RPI) modelling of wall boiling, and measured bubble departure diameters in conditions of the work in[20].

Figure 4 .
Figure 4. Comparison between measured wall temperatures, their values predicted using manual evaluation of the RPI model and boiling heat transfer correlations, for a typical boiling validation case.From[15].

Figure 4 .
Figure 4. Comparison between measured wall temperatures, their values predicted using manual evaluation of the RPI model and boiling heat transfer correlations, for a typical boiling validation case.From[15].

Figure 5 .
Figure5.Volume fraction distribution predicted with computational fluid dynamics (CFD) simulation using RPI modelling of wall-boiling.Note the mismatch between the size of the simulation domains, more than one metre in length, and the likely much smaller characteristic size of the steam bubbles generated at the heated wall, of perhaps only 10-100 micron in diameter.The computational mesh used to resolve the volume fraction distribution indicated in the figure consists of cells much larger than the single bubbles, which are replaced by a corresponding continuous distribution of the vapour phase.The image shows results for two different benchmark cases[22] ('bar2′ and 'bar4′ in the figure), using in in the CFD model the bubble departure diameter correlations of Tolubinskiy et al.[23] (panels (a) and (c) and Kocamustafaogullari et al.[21] (panels (b) and (d).For the same test case, different correlations for computing the bubble departure diameter return radically discrepant volume fraction distributions.From Colombo et al.[4].

Figure 5 .
Figure5.Volume fraction distribution predicted with computational fluid dynamics (CFD) simulation using RPI modelling of wall-boiling.Note the mismatch between the size of the simulation domains, more than one metre in length, and the likely much smaller characteristic size of the steam bubbles generated at the heated wall, of perhaps only 10-100 micron in diameter.The computational mesh used to resolve the volume fraction distribution indicated in the figure consists of cells much larger than the single bubbles, which are replaced by a corresponding continuous distribution of the vapour phase.The image shows results for two different benchmark cases[22] ('bar2' and 'bar4' in the figure), using in in the CFD model the bubble departure diameter correlations of Tolubinskiy et al.[23] (panels (a,c) and Kocamustafaogullari et al.[21] (panels (b,d).For the same test case, different correlations for computing the bubble departure diameter return radically discrepant volume fraction distributions.From Colombo et al.[4].
and do not usually cover the range of parameter space typical of industrially relevant flow conditions.

Figure 7 .
Figure 7.An example of the accuracy of current bubble departure models used to generate input data for RPI-based CFD simulation of boiling flows.The blue squares indicate application [32] of the energy based model of Ardron et al. [30].The black lines indicate results of application of the Klausner

Figure 7 .
Figure 7.An example of the accuracy of current bubble departure models used to generate input data for RPI-based CFD simulation of boiling flows.The blue squares indicate application[32] of the energy based model of Ardron et al.[30].The black lines indicate results of application of the Klausner[25] and Yun[34] models, and of the empirical fit to the Klausner model by Sugrue et al.[28].From[32].

Figure 9 .
Figure 9. Example of interface-capturing simulation of boiling at the scale of the single steam bubbles, from Tryggvason et al. [43].Computational meshes used for this kind of simulation are required to resolve the details of the vapour-liquid interface (grey surface in the figure), typically resulting in cell sizes of a few micrometres in realistic boiling conditions.

Figure 9 .
Figure 9. Example of interface-capturing simulation of boiling at the scale of the single steam bubbles, from Tryggvason et al. [43].Computational meshes used for this kind of simulation are required to resolve the details of the vapour-liquid interface (grey surface in the figure), typically resulting in cell sizes of a few micrometres in realistic boiling conditions.

Figure 10 .
Figure10.Microlayer formation and depletion beneath a bubble growing at a heated wall, a typical example of the importance of variables associated with the solid surface.From Hänsch et al.[67,73].

Figure 10 .
Figure10.Microlayer formation and depletion beneath a bubble growing at a heated wall, a typical example of the importance of variables associated with the solid surface.From Hänsch et al.[67,73].

Figure 11 .
Figure 11.Panels (a) and (b) indicate, respectively, the bubble shape and heat flux distribution at the solid surface reconstructed from measurement.Panel (c) shows comparison of heat flux distributions at the solid surface beneath a bubble from interface-capturing simulation (solid lines) and from the experiment (squares).The dashed lines represent predictions of the heat flux from the solid surface to the liquid obtained with the physical model used for RPI simulation of boiling at the component scale.From Giustini et al. [69,74].

Figure 11 .
Figure 11.Panels (a,b) indicate, respectively, the bubble shape and heat flux distribution at the solid surface reconstructed from measurement.Panel (c) shows comparison of heat flux distributions at the solid surface beneath a bubble from interface-capturing simulation (solid lines) and from the experiment (squares).The dashed lines represent predictions of the heat flux from the solid surface to the liquid obtained with the physical model used for RPI simulation of boiling at the component scale.From Giustini et al. [69,74].5.3.3.Towards CHF Prediction: Modelling the Collective Behaviour of a Small Population of Bubbles

Figure 12 .
Figure 12.Simulation of the collective behaviour of steam bubbles in pool boiling [56] at a heat flux of 300 / , tracking the evolution of a few bubbles for two seconds until the bubbles merge and

Table 1 .
Summative table of mass transfer models for boiling simulation.