Abstract
The objective of this paper is to give an overview of the capabilities of Eulerian bifluid approach to meet the needs of studies for nuclear safety regarding hydrogen risk, boiling crisis, and pipes and valves maintenance. The Eulerian bifluid approach has been implemented in a CFD code named NEPTUNE_CFD. NEPTUNE_CFD is a three-dimensional multifluid code developed especially for nuclear reactor applications by EDF, CEA, AREVA, and IRSN. The first set of models is dedicated to wall vapor condensation and spray modelling. Moreover, boiling crisis remains a major limiting phenomenon for the analysis of operation and safety of both nuclear reactors and conventional thermal power systems. The paper aims at presenting the generalization of the previous DNB model and its validation against 1500 validation cases. The modelling and the numerical simulation of cavitation phenomena are of relevant interest in many industrial applications, especially regarding pipes and valves maintenance where cavitating flows are responsible for harmful acoustics effects. In the last section, models are validated against experimental data of pressure profiles and void fraction visualisations obtained downstream of an orifice with the EPOCA facility (EDF R&D). Finally, a multifield approach is presented as an efficient tool to run all models together.
1. Spray Modelling and Wall Condensation in the NEPTUNE_CFD Code
A severe accident in a PWR nuclear power plant generally originates from a lack of cooling of the core, whose residual power can no longer be evacuated. In a few hours, due to multiple failures, human or hardware, including the failure of backup procedures, the structure of fuel elements deteriorates. Hydrogen is produced from the oxidation of zirconium sheaths and of structures of fuel elements during the phase of core degradation.
The hydrogen and water vapor thus produced are transferred to the containment and then transported by convection loops. Given the significant differences in density between hydrogen and other gases in the containment (nitrogen, oxygen, water vapor, carbon monoxide, carbon dioxide, etc.), hydrogen can accumulate preferentially in the upper parts of the compartments of reactor building. In case of strong heterogeneity, hydrogen can achieve high local concentrations that exceed the threshold flammable gas mixture.
Among the different safety systems for limiting the pressure increase during the course of the accident and the impact of possible combustion (deflagration), French and German PWR reactors use two types of mitigation means:(i)The passive autocatalytic recombiners (PAR): their role is to proactively oxidize hydrogen for preventing its accumulation in the containment. The catalytic recombiners initiate a controlled combustion, which is similar to a slow deflagration.(ii)Sprinkler systems: the injected water droplets cool the containment and lower the pressure by condensing steam on the droplets. They also promote mixing of gas by breaking quickly possible stratifications of the lightest gases.
The walls of the containment building and metal structures also play an important role from a thermal viewpoint. The walls, significantly cooler than the gas, condense the water vapor in the gas mixture and thus limit the pressure increase in the containment. Furthermore, the temperature difference between fluid and walls generates convection loops, enhancing the mixing of gases having different density.
In this first part dedicated to the hydrogen risk, our strategy is to increase gradually the level of complexity and to validate all the physical phenomena involved separately. The first step is the analysis and the understanding of gas stratification and transport phenomena. In order to achieve this, CFD calculations are performed with two kinds of solver, namely, NEPTUNE_CFD and Code_Saturne, for the PANDA test 17. NEPTUNE_CFD is dedicated to the simulation of incompressible and compressible multicomponent/multiphase flows, whereas Code_Saturne is dedicated to homogeneous incompressible or low Mach number compressible flows, with only one momentum equation representing the mixture of gases, liquid, and particles. NEPTUNE_CFD is mainly used for nuclear engineering, whereas Code_Saturne is used for nuclear and fossil energy engineering and for environment (as geophysical flows).
The vapor condensation along the wall of the containment building is the next step towards increasing complexity. The modelling proposed in NEPTUNE_CFD assumes that liquid droplets form along the wall. Vapor condensation on droplets makes them to grow. Once the droplet diameter reaches a critical value, gravitational forces compensate surface tension force and then droplets slide over the wall. Droplets can also join the surrounding droplets and form a film layer. Condensation terms implemented in Code_Saturne are calculated by correlations. Both approaches (single- and two-phase flow) are compared against the PANDA test 25 in this part. It is of relevant interest to be able to propose and validate several approaches. Indeed, this strategy allows improving the reliability of the results in the simulation of realistic scenarios, such as in TOSQAN ISP-47 [1].
The last step towards increasing complexity concerns the effects of spray modelling in the containment building. The two major thermal-hydraulic phenomena involved in spray behaviour in such applications are the thermodynamical effect of a spray (steam condensation on droplets and evaporation) and the dynamical effect (entrainment and mixing). But, steam condensation also leads to a local increase in the hydrogen concentration. Because of this local negative effect, dedicated CFD tools are required.
In the frame of the Severe Accident Research Network (SARNET), comparisons codes-experiments have been performed to evaluate the efficiency of these spray systems in terms of depressurization, hydrogen mixing, and radioactive aerosols scavenging for applications concerned by nuclear reactor accidents. We briefly recall the main steps in the following.
Firstly, tests on a single falling droplet under typical gas conditions representative of a nuclear accident have been proposed in the frame of the SARNET2 network [2, 3]. The drop diameter is modified by steam condensation or evaporation depending on the gas and droplet thermodynamical conditions.
Code-experiment and code-code comparisons are presented in [2, 3]. It is shown that the mass transfer terms are responsible for most of the differences observed between CFD codes and the kind of test that the errors can either compensate or cumulate. Since the errors are propagating over the droplet height fall, they could be nonnegligible for real containment cases. We will focus on simulations performed with NEPTUNE_CFD in this section.
Secondly, two tests have been produced in the TOSQAN facility in order to study the spray behaviour under severe accident conditions: TOSQAN 101 and TOSQAN 113 [2, 3]. The TOSQAN facility is a large enclosure (7 m3) devoted to simulate typical accidental thermal-hydraulic flow conditions in nuclear-pressurized water reactor (PWR) containment. Results obtained against these two cases have been carried out in the frame of the SARNET network.
The test is composed of a single nozzle producing a full-cone water spray [4]. The TOSQAN spray test 113 focuses on the dynamical effects of spray systems on entrainment and mixing of gases. The gas mixture is initially stratified: helium in the upper part of the vessel and air in the lower part of the vessel. As usually, helium is a simulant fluid used instead of hydrogen. When the spray is activated, gas mixture is entrained, generating a global flow that leads to a global mixing in the vessel after about 250 seconds. A faster mixing time was predicted by all the CFD codes.
As concerns test 101, an initial pressurization in the vessel is performed with superheated steam up to 2.5 bar. Then, steam injection is stopped and spraying starts simultaneously at a given water temperature (around 25°C) and water mass flow rate (around 30 g/s). It was found that some discrepancies were localized in the spray region which can thus have a high impact on the global results, since most of the heat and mass transfers occur inside this region.
Moreover, most of codes underestimated the transient pressure decrease even if the steady state pressure was correctly predicted. Let us examine this point that will be of relevant interest in the following sections. Figure 1 displays the number of moles of gas plotting against time [5–7]. The variation in the number of moles of gas in the vessel was equal to the variation in the number of moles of steam (the number of moles of air was constant). Let us consider the first phase beginning with spraying and lasting about 100 s. It was characterized by a sudden drop in the average gas temperature in the vessel. It can be seen that the number of moles of gas in the vessel rises as the gas temperature drops. It can be deduced from this that droplet evaporation was the main cause of gas cooling, with the droplets drawing from the gas the energy they required to evaporate. But, the air-vapor mixture was initially assumed as homogeneous in the benchmark exercise. With this assumption, thermal-hydraulics conditions of the test cannot lead to droplet evaporation. Because of mass and energy conservation, the steady state is finally correctly reproduced. We have to keep in mind that code-experiment comparisons may lead to discrepancies due to weaknesses in modelling but also because of exaggerated simplifications of initial and boundary conditions in calculations. Results obtained with NEPTUNE_CFD code are detailed in Mimouni (2010).
A step further in the understanding of the phenomena involved in an industrial geometry can be achieved in the numerical simulation of two interconnected vessels, which is a main objective of this section. The two experiments analysed with NEPTUNE_CFD in this section (ST3_0 and ST3_2) were conducted in the frame of the OECD/SETH2 project using the PANDA facility built at PSI (Switzerland) [8]. We will address in the following how these experiments of paramount importance provide valuable results into the formation of potentially dangerous hydrogen mixtures with increased hydrogen content in regions far away from the spray.
Moreover, benchmark calculations on a real PWR spray nozzle (CALIST) have been investigated with different codes, mainly CFD codes, from different institutions (Malet et al., 2015).
It is found that droplet velocities are almost well calculated 1 m below the spray nozzle, even if the spread of the spray is not recovered.
In this section, a two-fluid multidimensional simulation with NEPTUNE_CFD, on the basis of two interacting real PWR spray nozzles, is compared to the results obtained on the CALIST facility [9–12]. The influence of collisions between droplets is taken into account with a statistical approach based on the various outcomes of binary collisions.
1.1. The Numerical Solver and Physical Modelling
The solver, based on a pressure correction approach, is able to simulate multicomponent multiphase flows by solving a set of three balance equations for each field (fluid component and/or phase) [13]. These fields can represent many kinds of multiphase flows: distinct physical components (e.g., gas, liquid, and solid particles); thermodynamic phases of the same component (e.g., liquid water and its vapor); distinct physical components, some of which split into different groups (e.g., vapor and several groups of different diameter droplets); different forms of the same physical components (e.g., a continuous liquid field, a dispersed liquid field, a continuous vapor field, and a dispersed vapor field). The solver is implemented in the NEPTUNE software environment [5–7, 14–16], which is based on a finite volume discretization, together with a collocated arrangement for all variables. The data structure is totally face-based which allows the use of arbitrary shaped cells (tetrahedra, hexahedra, prisms, pyramids, etc.) including nonconforming meshes.
1.2. Governing Equations
The two-fluid model constitutes the following balance equations.
Mass Balance Equationswhere is the time, , , denote the volumetric fraction of phase , its averaged density, and velocity (-component), and is the interfacial mass transfer per unit volume and unit time. The phase index takes the values 1 for the continuous phase (gas) and 2 for the dispersed phase (droplets). In case of sectional approach as discussed in the last part, for the dispersed phase.
Momentum Balance Equationswhere is the pressure, is the gravity acceleration, is the interfacial momentum transfer per unit volume and unit time, and and denote the molecular and turbulent stress tensors (Reynolds stress tensor).
Total Enthalpy Balance Equationswhere is the phase-averaged enthalpy for phase and is the interfacial-averaged enthalpy. We assumed that the two phases are governed by the same averaged pressure field and we make no distinction between the pressures in the two phases or between the bulk pressure and the interface pressure for simplicity. The term denotes the interfacial transfer term of heat, the quantity being the interfacial area concentration.
The lift and added mass force can be neglected in our calculations. Thus, the interfacial transfer of momentum is assumed to be equal to the drag force. The wall friction terms do not appear in the momentum balance equations because solid walls are only present at the boundary of the flow domain and the wall friction is expressed through the wall boundary condition. The term denotes the wall-to-fluid heat transfer per unit volume and unit time for each phase. The two terms and denote the molecular and turbulent heat fluxes inside phase .
In the simulations described later, gas turbulence is associated with the -ε model, whereas dispersed phases turbulence is modelled with the Q2–Q12 model [17].
Equation of Transport on the Density of Drops. If is the density of drops (number of drops per unit volume), we write the following equation of transport of the density of drops:where is the source term related to fragmentation and the source term related to coalescence.
We neglect the coalescence and fragmentation phenomena in calculations of heat and mass transfer for a single droplet or for only one spray. Therefore, the variations of the droplet diameter are only due to the mass transfer between the droplets and the gaseous mixture (evaporation and condensation). Droplets are supposed to be of spherical shape:For interacting sprays, coalescence and fragmentation phenomena should not be neglected. Thus, in this case, a monodispersed approach may fail.
Mass Balance Equation of the Noncondensable Gas. The gas phase is a mixture of vapor and noncondensable gas. In case of air-vapor mixture gas, we solve the following equation:with and , the diffusion velocities of the air and vapor in the mixture, respectively, given by Fick’s law. The mass fraction of air and vapor is, respectively, and .
If another noncondensable gas is added to the gaseous mixture like helium, we solve in addition the corresponding mass balance equation with (mass fraction of helium) instead of .
Finally, the vapor mass fraction is calculated by considering that the sum of mass fractions is equal to 1.
Interfacial Mass Transfer Terms. The interfacial mass transfer can be summarized bywhere Sh and Nu are, respectively, the Sherwood number and the Nusselt number and are calculated by the correlations of Frössling/Ranz-Marshall [18]:where Sc, Pr, Re, and are, respectively, the Schmidt number, the Prandtl number, the Reynolds number, and the gas mixture diffusion coefficient.
is the film temperature and is equal to .
From (7) and (3), we can deduce that the interfacial transfer term of heat is negligible, excepted when the mass transfer term tends to zero. Thus the main phenomenon is the mass transfer: droplet evaporates or the vapor condensates onto the droplet surface. In the latter case (saturation ratio of steam in the vessel is about 1), pressure drops in the vessels because the number of moles of steam drops. The gas mixing follows an ideal gas law which implies that pressure and temperature are proportional.
As a consequence, the droplet evaporation is the main cause of gas cooling, with the droplets drawing from the gas the energy they required to evaporate water.
When tends to zero, gas and droplet temperature evolves under the action of , altering the equilibrium between and , which means that the mass transfer starts again.
In the case of droplet evaporation (low value of saturation ratio of steam in the vessel), a high rate of heat absorption is needed. Thus the gas temperature drops. By the ideal gas law, we deduce that pressure also drops.
Moreover, from (8), we conclude that heat and mass transfer depends on the Reynolds number and then on the relative velocity between gas and droplet. This dependence may cause discrepancies if it is not calculated accurately. Unfortunately, gas and droplet velocity are not available in most of experiments.
In case of vapor condensation at the wall, the modelling proposed in the paper assumes that liquid droplets form only heterogeneously at nucleation sites. Vapor condensation on droplets makes them grow. Once the droplet diameter reaches the critical value calculated, gravitational forces compensate surface tension force and then droplets slide over the wall. Droplets can also join the surrounding droplets and form a film layer.
From the formation of droplets by condensation, we deduce the mass transfer [14–16]: and so with being the surface density of drops on the wall. The number of droplets at the wall by surface unit iswhere is the number of droplets by surface unit formed at nucleation sites by vapor condensation and is the number of preexisting droplets by surface unit. We assume thatwhich is the formulation at the CFD scale of condensation properties experimentally observed at the droplet scale [19]!
1.3. PANDA Facility: Test 17
PANDA 17 is a test case which deals with the slow destruction of a stratification [20]. Initially Vessel-1 presents an Air-He mixture layer at ambient temperature and pressure (Figure 2). Then an air plume of 303 K is injected inside the facility to impact the layer and lead to its slow erosion. As the process is running, over the initial layer we can see a decreasing helium concentration whereas we can observe an increase of helium concentration on lower part and then a stabilization. The vessel temperature is constant and does not allow wall condensation. This experiment aims to test CFD code capacity to model correctly buoyancy effects. In this section, Code_Saturne and NEPTUNE_CFD codes are used in the calculations. Calculations obtained with Code_Saturne (named “single phase” in figures) in a previous work were in good agreement with experimental data [14–16]. Calculations performed with NEPTUNE_CFD are named “two phases” in figures.
In this test case, we compare the helium concentrations along time at different locations (Figure 3). Comparisons between experiment and calculations show a reasonable agreement qualitatively and quantitatively. We do not reproduce correctly the rude decrease of helium concentration on probe series 14 at elevation 5.626 m which is located close to the injector and the initial layer (Figure 4). This decrease is representative of the plume impact on the layer. This one influences the helium concentration evolution on the probes of the series 14. This is a consequence of our mesh which is coarse. On series 18 and 20, helium concentration evolution is well represented (Figures 5 and 6). The technical method to set up the case in two-phase flow approach is not based on mass sink/source terms. As a consequence and being associated with the fact that we use a coarse mesh, we have low frequencies oscillations on some figures compared to the single-phase computation on the same grid refinement. Further works should be done on mesh refinement and methodology to improve results on the impact of the plume. Nevertheless, those first results give us reasonable confidence for industrial purposes and give an idea of the degree of maturity of CFD codes for this kind of situation.
1.4. PANDA Facility: Test 25
The walls of the containment building and metal structures also play an important role from a thermal viewpoint. The walls, significantly cooler than the gas, condense the water vapor in the gas mixture and thus limit the pressure increase in the containment. Furthermore, the temperature difference between fluid and walls generates a convection loop, enhancing the mixing of gases of different density. The original approach presented in this work is the standard model in the code and is currently used in the simulation of the whole containment in industrial scenarios.
The test 25 experiment consists of injection of steam and helium through a vertical pipe of 0.2 m diameter in Vessel-1 (DW1) [21]. The mixture is vented into the PANDA Wetwells (WW1 and WW2) (Figure 7), which are used as a buffer volume smoothing the pressure increase in vessels DW1 and DW2 due to the gas injection. The conditions of the test 25 experiments are summarized in Table 1.
Initial Conditions. Pressure = 1.3 bar; fluid temperature = 27°C; structure temperature = 27°C; and gas composition contains air only.
The beginning of the transient is driven by positive buoyancy effects, for which the lightest fluid (helium) rises rapidly in the upper part of the two chambers. Then, as steam is condensed upon the cold structures surrounding the gases, the helium concentration in the first chamber relatively increases according to the other gases and the fluid becomes lighter than the injected fluid mainly composed of steam. A negative buoyancy effect is then observed, for which the injected fluid falls down in the first chamber. The sedimentation of the fluid in the first chamber is observed a short time after the beginning of the transient.
The vessel temperature evolution represented in Figure 9 is well reproduced but a little bit below the experiment results. Helium concentration evolution at higher and middle (Figures 12, 13, 14, and 15) locations is in good agreement with experiment. But at the bottom (Figures 10 and 11) of the two vessels, we obtain differences which can have different explanations such as residual phases or nonequilibrium between gravity and pressure gradient, for example. This observation should be linked with the vessel temperature evolution at the lower probe. As we have a high concentration of helium on the lower part of Vessel-1, wall condensation phenomena cannot occur due to the lack of vapor. Subsequently as the condensation flux is low, vessel temperature cannot increase.
Overall, results are reasonably correct and give an idea of the degree of maturity of CFD codes for this kind of situation.
1.5. HMT Benchmark
The IRSN CARAIDAS experimental setup [22] was used to study drop evolution under representative conditions of postaccident atmosphere. The cylindrical enclosure is of 5 m height and 0.6 m inner diameter. Homogeneous conditions are obtained with gas temperatures from 20 to 160°C, absolute pressures from 1 to 8 bar, and relative humidities HR from 3 up to 95%. Drop diameter optical measurements are performed at 3 elevations: at the top (drop generator, = 5 m), at = 2.49 m, and at the bottom ( = 0.61 m). The tests conditions are given in Table 2: “evaporation” and “condensation” tests are called EVAP- and COND- tests, respectively. Let us precise that only one drop is ejected in each case and is characterized by its initial temperature () and its initial vertical velocity (). Benchmark specifications are given in Malet (2013). More details of the experiments can be found in [23].
Figures 16 and 17 show that the mean relative error is less than 5% and we can reasonably conclude that calculations are in good agreement with experimental data.
In case EVAP24, heat and mass transfer are cancelled when the droplet diameter decreases under a specified value in order to avoid dividing by zero when the droplet evaporates. Results proved to be sensitive to the minimum value of droplet diameter prescribed in this case (Figure 17). This numerical problem should be studied in further calculations.
Mesh and time convergence have been carefully investigated (Figure 16): 15 cells in radial direction and 125 cells in the axial direction for the coarse mesh; 30 cells in radial direction and 250 cells in the axial direction for the medium grid; and 60 cells in radial direction and 500 cells in the axial direction for the fine grid. Nevertheless, dealing with single droplet in an eulerian approach proved to be a challenging task regarding the boundaries conditions on liquid volume fraction at inlet. As a consequence, these cases illustrate the well-known problem of vanishing phases which remains an important topic in mathematical journals.
1.6. PANDA Facility: ST30 and ST32
The goal of the tests is to study the effect of spray on the containment depressurization and on the distribution of noncondensable gases for a generic, but still containment relevant configuration. Each test features initially a helium containing layer at the top of Vessel-1. The two vessels, each has a height of about 8 m and a diameter of 4 m, are connected with a 1 m diameter interconnecting pipe (IP).
The spray nozzle is oriented vertically downward in Vessel-1 and its outlet is located at 6.9 m from the bottom of Vessel-1. The nozzle has 6.4 mm outlet diameter and produces a conical spray pattern with an opening angle of 30° for the specified flow rate and pressure.
Two containment spray tests have been performed. Denoted as ST3 0 and ST3 2, these experiments differ only in the composition of the initial gas atmosphere. For ST3 2, the helium-rich layer was created in an initial steam-air environment at the top of Vessel-1, and the remaining volume of Vessel-1 as well as the full volume of Vessel-2 was filled with the steam-air mixture. This test is conducted at hot conditions and at an elevated pressure.
The cold test ST3_0 was performed with air and helium in order to exclude heat transfer between the spray droplets and the ambient gas and steam condensation. It serves as a reference test where the gas mixing and break-up of the initial layer are only driven by the hydrodynamic interaction between droplets and ambient gas. For both cases, the gas mixture is initially stratified: helium in the upper part of Vessel-1 and air or air-vapor in the lower part of Vessel-1. For both cases, the water injection flow rate is 840 g/s. The water injection temperature is 40°C for ST3_2 and 25°C for ST3_0.
Boundary and initial conditions are presented in Figures 18–20. Unlike TOSQAN 101 test, molar fractions and gas temperature are not assumed to be homogenous initially (Figures 19-20).
For ST3_2, the initial pressure is a typical postaccident containment vessel pressure and the initial gas temperature corresponds to the saturated vapor temperature at this pressure (Figure 18).
Details about PANDA facility and instrumentation can be found in Erkan et al. [8].
We neglect the collision between droplets in the calculations because only one spray is activated. A monodispersed approach is chosen and the Sauter mean diameter at injection is taken equal to 0.582 mm. The meshing contains 22000 cells.
Figure 21 shows the helium molar fraction after 171 s and 4000 s for ST3_0. At 171 s after the spray activation, the initial helium stratification can be observed in the upper part of Vessel-1. Moreover, a stratified mixture flows through the IP towards Vessel-2. Since an air-helium mixture with higher helium content and the same temperature is lighter than pure air, the air-helium mixture coming from Vessel-1 travels along the top of the IP in accordance with experimental measurements [8]. At = 4000 s, the gas mixture is uniformly mixed in Vessel-1 whereas a helium-rich mixture fills Vessel-2 from the top.
The time-dependent helium molar fractions in Vessel-1 and Vessel-2 are displayed in Figures 22 and 23. The position of the probes can be deduced from Table 3 and Figure 8. The good agreement with experimental data shows that the dynamical effects of spray systems on entrainment and mixing of gases are correctly taken into account. But a faster mixing time was predicted in TOSQAN 113 test. According to the authors, the gas entrainment created by the spray interacts in a lesser extent with the flow generated at wall in PANDA geometry. Thus, the gas flow is simplified in the ST3_0 test and is easier to simulate.
(a)
(b)
Without heat and mass transfer, the momentum exchange between droplet and gas is correctly predicted, which is of relevant interest to evaluate discrepancies for ST3_2 test where mass, momentum, and energy exchanges occur.
Figure 24 shows the helium molar fraction after 171 s, 1001 s, and 4000 s for ST3_2. As in the previous case, firstly, an air-steam-helium mixture coming from Vessel-1 travels along the top of the IP towards Vessel-2. But, after the spray activation, vapor condensates onto surface droplet and thus vapor concentration decreases in Vessel-1, as described in Figure 25. As a consequence, the air concentration increases in the gas mixture whereas the gas temperature drops in Vessel-1 as explained in section 0 and described in Figure 26. At the beginning, to equilibrate the mass between Vessel-1 and Vessel-2, an air-steam mixture coming from Vessel-2 travels along the bottom of the IP towards Vessel-1 and then reinforces the air concentration in Vessel-1.
The gas mixture becomes heavier and colder in Vessel-1 than the gas mixture in Vessel-2. As a result, the helium-rich mixture coming from Vessel-1 travels along the bottom of the IP (Figure 24 at = 1001 s) in accordance with experimental measurements [8]. Figure 24 at = 4000 s shows that the helium-rich mixture travelling the IP fills Vessel-2 from the bottom, which was not obvious and makes this experiment of relevant interest. During this phase, a hot air-vapor rich mixture coming from Vessel-2 travels at the top of the IP towards Vessel-1.
An accurate description of the physical phenomena involved during the transient of tests ST3_0 and ST3_2 can be found in [8]. Regarding test ST3_2, main phenomena are in accordance with experimental data but the time-dependent pressure in Vessel-1 (Figures 27 and 28) shows that an important phenomenon is not modelled which makes results unacceptable.
By considering the time-dependent gas mixture density in Vessel-1, Erkan et al. [8] note that helium layer erosion develops in tests ST3_2 slightly later (see Figures 21 and 24 at = 171 s). Indeed, Figures 22 and 29 show clearly the delay and make calculations and experimental data in good agreement. But, we disagree with Erkan et al. [8] and we consider that spray is inefficient on the helium layer erosion because at the beginning the spray does not exist practically: experimental profiles suggest that droplets evaporate.
Erkan et al. [8] have measured the time-dependent saturation ratio of steam in Vessel-1. The initial values are relatively low but increase throughout all measurement locations after the spray activation. This latter experimental observation suggests also that droplets evaporate.
By considering the schematic view of PANDA facility in Figure 18, we see that spray is injected at 1.2 m from the top of Vessel-1, initially maintained to 129°C. Thus, we can suppose that the injection pipe is at a temperature between 129°C and 40°C (water temperature injection). To test the assumption, we have performed a simulation with a water temperature injection equal to 90°C for 100 s and 40°C after that. Results are presented in Figure 28: the experimental data support the assumption of hot water injection at the beginning of spray activation. As a consequence, with this assumption, time-dependent vapor profiles are really in better agreement with experimental data (Figure 30).
The time-dependent helium molar fractions in Vessel-1 are displayed in Figure 29. The good agreement with experimental data is encouraging in the frame of hydrogen risk for containments applications.
Finally, heat, mass, and momentum transfer proved to be correctly predicted for only one spray.
For several sprays, can we neglect the coalescence and fragmentation phenomena? Furthermore, at outlet of the nozzle, can we continue to describe the experimental size distribution of droplet diameter by a monodispersed approach? These issues are addressed in the following section.
1.7. Modelling and Validation of Droplet Polydispersion and Collisions
The objective of this section is to propose a generalization of the monodispersed approach described and validated in the previous sections.
1.7.1. Experimental Measurement of PWR Containment Spray Characteristics
Experiments have been carried out at the French Institute of Radiological Protection and Nuclear Safety (IRSN) and half were financially supported by Électricité de France (EDF R&D), on the CALIST facility (Characterization and Application of Large and Industrial Spray Transfer [24]). In a room of 7 × 6 × 3.5 m3 dimensions, the setup is composed of a supplying hydraulic circuit and, for these experiments, of two interacting spray nozzles with a flow rate of 1 l/s at a relative pressure of 350 kPa for each nozzle and separated by 42 cm. The water spray, with a temperature of around 15°C, is collected in a 5 m3 pool.
1.7.2. Characteristics of Droplets at 20 cm from the Nozzle
Measurements performed at 20 cm from the nozzles are used as inlet conditions of the numerical simulations. These nozzles are used at a relative pressure of 3.5 bar, for a mass flow rate of 1 kg/s. At this distance, due to the hollow cone created by these nozzles, most of the droplets are concentrated in an annular area located between 8 cm and 15 cm from the nozzle axis (), with a maximum of presence at 11 cm.
The geometric mean diameter (), Sauter mean diameter (), and mean velocities are displayed in Figure 31 as functions of the distance from the nozzle axis. varies between approximately 240 μm and 330 μm. varies between 360 μm and 520 μm. This implies dispersion in size. The axial velocity is maximum close to the nozzle axis: it is 20 m/s at 8 cm and then decreases radially to 13 m/s at 15 cm. The radial velocity is maximal far from the nozzle axis and equals 7.7 m/s. The azimuthal velocity is very low and varies between 0.17 m/s and 0.34 m/s. This means that the swirl created by the nozzle is attenuated very quickly in the first centimetres when atomization occurs.
Figure 32 shows the local spray size and axial velocity distributions. It can be noticed that the shape of the size distribution does not depend on the distance from the nozzle axis. The size distribution can be approximated with a log-normal law [24]. The repeatability is very good for the and axial velocity measurements. Uncertainties are higher for radial and azimuthal velocities, because direct measurements of these two values are not possible with our PDI.
1.7.3. Modelling of Droplet Polydispersion
Greenberg et al. [25] developed a method to model particle polydispersion in size in Eulerian simulations. The idea was to consider the dispersed phase as a set of continuous “fluid” media: each “fluid” corresponding to a statistical average between two fixed droplet sizes, namely, a section. The spray was then described by a set of conservation equations for each “fluid.” In our case, both interacting sprays will be considered as independent fluids, and for each spray or fluid, size distributions are divided into sections of fixed diameter. Sections are chosen as fixed in size, and they exchange mass and momentum in order to model evaporation/condensation or collision phenomena. Mass and momentum balance are solved for each section.
1.7.4. Modelling of Droplet Collisions
The several-fluids model constitutes a mass balance equation, where represents each section. The mass transfer in mass balance equation becomes , where is the mass transfer per unit volume and unit time due to collisions. It is assumed that no evaporation or condensation occurs in this section: .
constitutes a source term and a sink term :The following can be written:where is the mass transfer from sections to after a collision between droplet of class and droplet of class .
A supplementary term is added at the right hand side of the momentum balance equation and is given by , where is the velocity of the section resulting from the collision between and .
and can be calculated with a collision frequency and a modelling of collision issue. For this latter modelling, five binary collision regimes can be pointed out: bouncing, coalescence, reflexive separation, stretching separation, and splashing [26]. Looking at the collision pictures (Foissac, 2010), it is possible to determine the final daughter diameter as a function of the initial “parent” diameters, using mass conservation. These values are summarized in Table 4. For the splashing regime, a value of 20 droplets has been estimated, but it should be considered as a first approximation. All these collision issues can be represented by a graph depending on the Weber number and the impact parameter.
As shown in Rabe (2010), it is appropriate to define a symmetric Weber number starting from first mechanical principles. Using the momentum balance, and assuming that the droplets are spherical and have the same density, the symmetric Weber number is then expressed byRabe et al. [27] proposed simple formulae expressing the boundaries of collision outcomes fields as a function of the symmetric Weber number. The final equation for the critical impact parameter ( represents the impact parameter for which the transition between two regimes is observed) at which the transition between reflexive separation and coalescence occurs is thenBased on another balance of energies, the critical impact parameter for the transition between coalescence and stretching separation can be expressed asFor larger droplets and velocities, only separation regimes can be observed, namely, reflexion and stretching. The ratio of the reflexive kinetic energy and the stretching kinetic energy can then be written and the critical impact parameter is then derived:with , a viscous dissipation coefficient, found experimentally to be equal to 0.92 and with dimensionless number that is found to be 0.28 according to experimental results (Rabe, 2010). These three models describing the transition curves between collision regimes are described more in detail in Rabe (2010). They are valid under ambient gas conditions for droplet sizes between 200 and 400 μm, with velocities up to 10 m·s−1.
Based on an energy balance, Estrade et al. [28] proposed an equation for the transition to bouncing (where is the fraction of volume interaction and the diameter ratio):according to experimental results.
Finally, splashing is assumed to occur when symmetrical Weber number is higher than 20 which is the very first modelling that needs to be confirmed by experiments.
It is also necessary to evaluate the collision frequency between droplets from sections and . Pigeonneau and Feuillebois [29] proposed the following expression:where and , and , and , and and are, respectively, the diameter, the number concentration, the droplet kinetic energy, and the fluid-droplet velocity correlation coefficient of sections and . is the radial distribution function introduced by Patino-Palacios and Simonin [30] and is equal to 1 in the calculations.
Therefore, and are obtained from the product of the collision frequency between sections and into the probability of a collision outcome derived from Rabe et al. [27], into the mass or velocity difference between sections and associated with the collision outcome.
The sectional method and the polydispersion modelling have been validated against Direct Numerical Simulations (DNS) of particle clouds in homogeneous isotropic turbulence, without gravity, in a previous work [24]. The next section deals with the validation of the drift part of the collision frequency, since this is the main phenomenon responsible for the collision in the top of the reactor containment. Calculations have been carried out at EDF R&D and half were financially supported by IRSN.
1.7.5. Numerical Simulation of Two PWR Interacting Sprays
Interacting sprays, characterized on the CALIST facility, are simulated inside a parallelepiped mesh of 800,000 hexahedra regular cells, representing a domain of 1.20 × 0.80 × 2 m. All boundaries are considered as free outputs, except the top face which contains the input and walls around, with no friction, and where velocity can only be tangential. Since the spray produced by these nozzles is a hollow cone one to the location 20 cm from the nozzle, the input domain is modelled by two annular rings of 18 cm internal diameter and 26 cm external diameter. Droplets are injected from this annular ring with the size distribution presented in Figure 32. On this figure, it can be seen that, in function of the distance to the nozzle axis, the axial velocity decreases from 20 m/s to 15 m/s.
In these simulations, it was assumed that the injection velocity is independent of the position; it was chosen with a value of 18.6 m/s, that is to say the velocity at 11 cm from the nozzle axis, where the volumetric fraction is maximal. Estimating the value of the radial velocity for the simulation is more difficult. Indeed, this value is very important since it is the main component of the relative velocity of the droplets when spray interacts, and so is the value of the Weber number and the collision frequency. A value of 7.7 m/s was chosen according to the results presented in Figure 31. The azimuthal velocity was neglected due to its low value.
Each spray size distribution was separated in 9 sections (Figure 33), whose void fractions were adjusted from the assumed droplet size distribution so as to obtain a mass flow rate of 1 kg/s, as measured on the real PWR nozzle for a relative pressure of 3.5 bar.
Before simulating two PWR interacting sprays we compare monodispersion and polydispersion for one PWR spray. As previous cases, time and grid convergence have been carefully studied (Figure 34) which is of paramount importance. Without such verification, calculations of industrial configurations are subject to criticism.
Figures 36 and 37 show the axial and radial droplet velocity calculated with a monodispersed approach whereas Figures 38 and 39 deal with the polydispersed approach. Measurements performed at 20 cm from the nozzles are used as inlet conditions of the numerical simulations and are correctly taken into account in the polydispersed approach.
Inlet conditions are coarsely imposed in monodispersion which leads to discrepancies. Thus, polydispersion improves the prediction of droplet velocity and droplet diameter (Figure 35).
Experimental and numerical local size distributions obtained for two interacting sprays are compared in Figure 40 for different positions along the symmetrical axis. It is clear that the droplet size decreases since the mean geometric diameter is about 300 μm before spray interaction and about 200 μm after spray interaction (Figure 40). This decrease can have two origins. First, it can be due to collisions at high Weber number that occur when sprays interact: in the interaction area, collision frequency reaches a maximum of about 1011 collisions·m−3·s−1, and the Weber number is very high, so that collisions could lead to break-up.
This size decrease is also due to the entrainment of the smallest droplets in the direction of the symmetrical axis [31]. The smallest droplets are drifted away in the air flow (Figure 42), whereas the biggest droplets (Figure 41), having more inertia, are not altered in the spray interacting area. At this stage, we still have to separate the effects of these two phenomena.
Many parameters have to be tested in order to evaluate their influence. The first one is the radial velocity at the inlet, since it is involved in many critical parameters like the Weber number and the collision frequency. The difficulty is that its value is bound to an uncertainty in the measurement. The sensitivity to the choice of intervals of the size distribution is also part of the future work.
Figure 43 represents the mean Sauter diameter in several cases: with or without collision and with or without polydispersion. As already discussed, the monodispersed case leads to large discrepancies. The polydispersed cases give a reasonable agreement with experimental data. But, the splashing phenomenon should be improved in further calculations because it seems that too much small droplets are created.
1.8. Conclusion about Our Current Capabilities to Simulate Sprays
Gas stratification and vapor condensation at the walls in the containment building have been assessed with two kinds of CFD approaches. Calculations compare reasonably well with experimental data and give an idea of the degree of maturity of CFD codes.
The reliability and applicability of the spray model implemented in the NEPTUNE_CFD code have been also addressed in this section: heat and mass transfer (HMT) on single droplet (CARAIDAS experiment); spray behaviour under severe accident conditions in a closed cylindrical vessel with a volume of 7 m3 (TOSQAN experiment).
A step further in the understanding of the phenomena involved in an industrial geometry was achieved in the numerical simulation of two interconnected vessels using the PANDA facility built at PSI (Switzerland). Two experiments were analysed with NEPTUNE_CFD in this section: ST3_0 and ST3_2 test. These experiments have provided valuable results for the formation of potentially dangerous hydrogen mixtures with increased hydrogen content in regions far away from the spray.
Calculations were in good accordance with experimental values for ST3_0 test, where the main phenomenon was the momentum exchange between droplet and gas.
A qualitative agreement was achieved for ST3_2 in a first attempt. But measurements and considerations about the PANDA facility have suggested that during about 100 to 200 s the water temperature injection was hotter than 40°C which finally led to the droplet evaporation instead of vapor condensation onto droplets surface. Under this latter assumption, calculations were found in reasonable agreement with experimental data.
Finally, a numerical simulation of the interaction between two PWR containment sprays has been performed with a new model of polydispersion and collision of droplets. The droplet size and velocity distributions at a distance of 20 cm below the spray nozzle outlet have been precisely measured and used as input data in the calculations. The water droplet polydispersion in size has been treated with a sectional approach. The influence of collisions between droplets is taken into account with a statistical approach based on the various outcomes of binary collisions. A two-fluid multidimensional simulation, on the basis of two interacting real PWR spray nozzles, is compared to the results obtained on the CALIST facility and shows a good agreement. However, monodispersed and polydispersed approach gives similar results regarding the air entrainment caused by droplets. Furthermore, polydispersed approach is needed for refined mesh but, in containment applications, computational cells remain too coarse (several centimetres) because of calculations of long transients (several hours).
These results allow us to continue on sensitivity studies in order to evaluate the most important phenomena involved in the droplet characteristics evolution (condensation, evaporation, entrainment, and collision). The knowledge of these characteristics could be important to evaluate the efficiency of these spray systems in terms of depressurization, hydrogen mixing, and radioactive aerosols scavenging for applications concerned by nuclear reactor accidents.
2. Recent Advances in Nucleate Boiling Modelling and Application to DNB
When a liquid is flowing onto a heated wall, the heat transferred by the wall to the liquid leads part of the liquid to evaporate, therefore producing a two-phase bubbly flow along the wall. This kind of situation, called the nucleate boiling regime, is an efficient manner to evacuate the heat produced into the wall, for example, by Joule effect or by a nuclear reaction. In nucleate boiling, the heat flux increases and reaches a maximum value with increasing surface temperature. Unfortunately, further increase in the surface temperature results in decreasing heat flux as the transition from nucleate boiling to film boiling takes place. The maximum heat flux that can be obtained by nucleate boiling is referred to as the critical heat flux (CHF). In the case of controlled heat flux, a slight increase of heat flux beyond the CHF can cause the surface temperature to rise to a value exceeding the surface material’s maximum allowable temperature. This in turn can cause severe damage or meltdown of the surface.
As a consequence, CHF has been extensively studied in the last five decades, as a major limiting phenomenon for nuclear power plant capabilities, as well as in other industries.
The scope of the present work is limited mostly to the boiling crisis in subcooled-flow boiling, which is of interest in the design of fuel assemblies used in nuclear fission, pressurized water reactors (PWR). This kind of situation is named departure from nucleate boiling (DNB).
To predict CHF, many empirical correlations have been developed as well as a few theoretical models. Although empirical correlations can be very reliable in the range of conditions where they have been established, their use outside this domain is very hazardous. On the contrary, theoretical models, by taking into account the basic mechanisms involved in the CHF phenomenon, should better adapt to any new flow boiling configuration.
In a previous work, a new mechanistic model has been implemented in a computational multifluid dynamics tool leading to wall temperature excursion and onset of boiling crisis. Even if this model covers a large physics scope in terms of mass flux, pressure, quality, and channel diameter, it was necessary to improve the calculation of the critical heat flux in some cases as for high values of quality and low pressure at outlet.
As a consequence, this section aims at presenting the generalization of the previous DNB model and its validation against 1500 validations cases. Another way of improvement comes from the fact that most of nucleate boiling models implemented in the CFD codes are based on the calculation of the site density where bubbles are emitted. This quantity is evaluated by an empirical correlation which depends on thermal-hydraulics conditions and material surface. Such formulation leads to large uncertainties because of the small amount of data available on site density and the difficulty to measure it.
Furthermore, the wall temperature proved to be overestimated with the standard model.
Thus the second objective of this section is to propose a new mechanistic model devoted to CFD scale and allowing similar results for the thermal-hydraulics variables (velocity, liquid temperature, and void fraction near the wall) but improved results for the wall temperature.
Firstly, the general model we use for two-phase boiling flow simulations is presented. The next paragraph is dedicated to validation of DNB tests in a tube. A new nucleate boiling model is proposed. Finally, conclusions are drawn about our current capabilities to simulate DNB and perspectives for future work are given.
2.1. Governing Equations
In this section, we solve the same balance equations as already described in the previous section for both liquid and vapor phases.
The forces exerted on bubbles are the averaged drag, added mass, and lift and turbulent dispersion forces. The liquid turbulence is modelled by the -ε SSG model adapted to bubbly flows (Mimouni, 2010, 2011).
The bubble size distribution modelling has been developed for bubbly flow based on the moment density method [32], where we assume that, in a computational cell, all the bubbles have the same velocity and the same temperature despite possibly different diameters.
2.2. Wall Function for Boiling Flow
In subcooled-flow boiling, the liquid velocity profile in the boundary layer is significantly disturbed by the bubble formation and detachment mechanisms on the heated wall. In the literature, an overprediction of liquid and gas velocity distributions in the boiling boundary region has been reported. The use of single-phase wall law may be one of the main reasons for these results.
In order to take into account the influence of bubbles in the near wall area, a modified logarithmic law of the wall was suggested by Mimouni et al. [5–7] which is usually used for turbulent flows over rough walls and reads aswhere , , and ( is the wall shear stress). Here is the known velocity tangential to the wall and is the distance from the wall. Coefficients and are standard single-phase constants with the values of 0.41 and 5.3, respectively. The last term represents the offset of due to the wall roughness:where for sand-grain roughness and is the so-called roughness Reynolds number:where and , where α and are, respectively, the void fraction and the bubble diameter in the adjacent cell at the wall.
When the void fraction tends to zero, the wall law tends to the single-phase formulation. Furthermore, this relation depends on bubble diameter and bubble density at the wall (void fraction at the wall), which is physically expected. This formulation proved to be a key point in the CHF simulation and is the standard option in the use of the NEPTUNE_CFD code.
2.3. Standard Wall Transfer Model for Nucleate Boiling
The CHF phenomenon has received considerable attention in the past, and different mechanisms have been proposed to interpret its cause.
Weisman and Pei [33] assume that DNB appears when a thin layer adjacent to the wall reaches a limiting bubble concentration due to the inability of the main stream to remove the bubbles, because the turbulent eddies are too small to influence the bubble trajectories. The critical void fraction in the layer corresponds to a maximum packing of bubbles and is estimated at 82%.
Lee and Mudawwar [34] postulate that small vapor blankets are formed due to the piling of bubbles flowing along the wall after their departure. A dry spot may appear when the vapor blanket length is such that a Helmholtz instability of the liquid-vapor interface occurs. DNB is then assumed to appear when the vaporization rate overcomes the liquid flow in the sublayer.
In the present section, in a first, simplified approach, and following the analysis of Kurul and Podowski [35], the boiling heat flux is split into three terms:(i)A single-phase flow convective heat flux at the fraction of the wall area unaffected by the presence of bubbles(ii)A quenching heat flux from bubbles departing from the wall and bringing cold water in contact with the wall periodically(iii)A vaporisation heat flux needed to generate the vapor phase.
In order to determine the different terms of the wall heat flux, three quantities are required: the lift-off diameter, the nucleation frequency , and the nucleation site density Na. They are computed thanks to a set of correlations given in [35, 36].
This basic wall heat flux partitioning model assumes that the amount of water on the wall is sufficient to remove heat from the wall and to be used for evaporation. Superheating of the vapor that occurs at high void fractions is not modelled. Given all this, the basic heat flux partitioning model cannot be used under critical heat flux conditions. In order to take into account the phenomenon of temperature excursion at DNB conditions, the heat flux partitioning model has been generalized by adding a fourth part of the wall heat flux, , diffusive heat flux used to superheat the gas phase [37] (Mimouni, 2010):where is the wall heat transfer coefficient calculated from the temperature wall function for the vapor phase, is the vapor temperature at the centre of the wall-adjacent cell, and is the wall temperature. Thus, the heat flux imposed at the wall is written aswhere is a phenomenological function, which depends on the liquid volume fraction and takes care of the numerically smooth transition between the nucleate boiling regime and the boiling crisis regime such as .
This function is represented in Figure 44.
To ensure grid independence, the liquid temperature in the nucleate boiling model is calculated from the logarithmic temperature profile at the given nondimensional distance from the wall rather than from the centre of the wall-adjacent cell. is the nondimensional distance based on the lift-off bubble diameter calculated at the onset nucleate boiling and remained constant along the simulation. The same process is applied to the liquid velocity [37].
Moreover, if the liquid temperature in the nearest cell at the wall tends to the saturation temperature, then the single-phase-flow convective heat flux tends to 0, and the corresponding energy contributes to the vaporisation heat flux.
We can now define the DNB model [38]: If reaches the saturation temperature, then ; otherwise .
The condition can be seen as the condition to detect the interface between the liquid and the vapor film formed at the wall when the boiling crisis happens. The condition can be seen as a “Weisman criterion” [38].
2.4. Calculations of DNB Tests in a Tube
Among various tests, we have selected mostly negative-quality (at the outlet) tests such that the condensation effects dominate coalescence/break-up effects, the modelling of which is far from being reliable today. As a consequence, the bubble diameter drops with distance from the wall.
2.4.1. Calculation Procedure
The Russian Academy of Sciences produced a series of standard tables of CHF as function of the bulk mean water condition and for various pressures and mass velocities for a fixed tube diameter of 8 mm [39] as presented in Table 5.
For tube diameters other than 8 mm, the CHF is given by the approximate relationshipfor between 4 and 16 mm. Let us note that this approximate relationship induces a supplementary relative error as discussed in [38].
The flow is assumed to be axisymmetric; therefore two-dimensional axisymmetric meshing is used. The calculation is started with the wall heat flux equal to 70% CHF. The wall heat flux is then increased by 5% progressively. After each step, the wall heat flux reaches a plateau in order to stabilize the boiling flow. This procedure is repeated until the wall heat flux is equal to 130% CHF. In the calculations, CHF is detected when the wall temperature exceeds 1000 K. In fact, the zircaloy clad tubes start to degrade at about 1000 K. Because of the sudden rise in temperature, results are weakly sensitive to the wall temperature chosen for CHF detection.
2.4.2. Results and Generalization of the DNB Model
The sensitivity of mesh refinement has been evaluated in previous works [38] and was found acceptable.
In contrast to empirical relations, theoretical models implemented in the CFD code that take into account the basic mechanisms involved in the CHF phenomenon should better adapt to any new flow boiling configuration and to geometry. Thus, the DNB model has been assessed in a large scope in terms of mass flux, pressure, quality, and channel diameter in [38] and some results are reproduced in Figures 45 and 46.
As already discussed, we propose in this section an improvement of the DNB model. Indeed, regarding strongly subcooled cases with high values of liquid mass flowrate at inlet, it is quite difficult to reach the condition . The new DNB model can be simply seen as a generalization of the former one.
We keep the following condition (“volume condition”): If reaches the saturation temperature, then ; otherwise ,
which is needed in the definition of .
And we add a supplementary condition (“surface condition”): If reaches the saturation temperature, then ; otherwise , where is the wall area affected by bubbles nucleation. This condition leads to a new function .
Finally, the heat flux imposed at the wall is written aswith .
This model is validated against 1500 validation cases: pressure between 100 and 160 bar, quality between −0.3 and 0.3, tube diameter of 6, 7, 8, 9, and 10 mm, and liquid mass flow rate between 2000 and 5000 kg/m2 s. Figures 47, 48, and 49 present some results as a function of the quality . Figure 50 summarizes the results and shows a reasonable agreement with the experimental data.
Linear alignment of experimental and calculated CHF in Figure 50 suggests nonstochastic (bias-type) errors in calculations. Discrepancies seem to increase for high CHF values. Indeed, the wall heat flux is increased by 5% progressively and maintained constant until a stationary state is reached. If the CHF = 1 MW/m2, the wall heat flux step is equal to 0.05 MW/m2 but equal to 0.5 MW/m2 for CHF = 10 MW/m2 which leads to a coarse discretization regarding the wall heat flux. As a consequence, high CHF values are less well predicted.
2.5. New Nucleation Boiling Model Devoted to High Pressure Flows
We have improved the DNB model in the previous section by taking into account not only the void fraction in the cell but also the wall area affected by bubbles. Let us note that results obtained with the new DNB model are similar and sometimes slightly degraded compared to the former model. But the new DNB model covers a larger physics scope in terms of mass flux, pressure, and quality.
Another way of improvement comes from the fact that most of nucleate boiling models implemented in the CFD codes are based on the calculation of the site density where bubbles are emitted. This quantity is evaluated by an empirical correlation which depends on thermal-hydraulics conditions and material surface. Such formulation leads to large uncertainties because of the small amount of data available on site density and the difficulty to measure it.
Furthermore, the wall temperature proved to be overestimated with the standard nucleation boiling model. Thus the objective of this section is to propose a new nucleation boiling model devoted to CFD scale and allowing similar results for the thermal-hydraulics variables (velocity, liquid temperature, and void fraction near the wall) but improved results for the wall temperature.
2.5.1. Presentation of the Model
In the previous model, the code was using a loop on the heat flux equation to find the wall temperature by an iterative process. The novelty we proposed is to calculate the wall temperature as a function of the imposed wall heat flux; then an iterative process is carried out on the nucleation site density to solve the heat flux balance. The key idea is that the piece of experimental data used to determine a correlation for the wall temperature is much more important than that for nucleation site density. Furthermore, temperature measurements are more reliable than those of the nucleation site density. These facts lead to a correlation more general and more robust for wall temperature.
The range of application covered by the Liu-Winterton correlation [40] is very large for pressure, heat flux, and mass flow and as a consequence has been selected in the following:where(i) is forced convective heat transfer coefficient given by Dittus-Boetler [41],(ii) is nucleate boiling heat transfer coefficient given by Cooper [42],(iii) is factor to amplify the heat transfer by forced convection because of the turbulence induced by the bubble,(iv) is reduction factor of the heat transfer by nucleate boiling because of the convective flow:where Pr is Prandtl number and Rel is Reynolds number for the liquid phase.
The Liu-Winterton formula has been adapted to CFD approach by considering the liquid temperature at the nearest cell at the wall. Moreover, for each adjacent cell at the wall, we determine the nearest cell adjacent at the opposite wall. is calculated by using the maximum value of the liquid velocity along the segment [] by using a diffusion process.
First of all, bubble detachment diameter and bubble detachment frequency are determined outside the loop. Then, the wall fraction occupied by bubble nucleation is calculated by and next the heat fluxes.
For the boiling heat flux, the flux balance is used:Then, by using the usual expression of the boiling heat flux [35], the nucleation site density comes as follows:where is the latent heat and is the bubble detachment volume.
Thanks to the new density, the new influenced area and heat fluxes are calculated. If the influenced area is converging and if the calculated boiling heat flux is converging, the calculation goes on; if not, it keeps iterating.
Let us note that this approach does not remove all the physics-based approach of the Kurul-Podowski heat partitioning model which is still valid and used in the calculations. Let us underline that, thanks to an iterative process, only the site density correlation is just replaced by a more general and more robust correlation.
2.5.2. Validation Cases
Boiling Flow in a Vertical Pipe, with Water as Working Fluid. This Water-Steam test case consists in the stimulation of a nonequilibrium boiling flow as created by Bartolomei et al. [43]. In the three experiments, the outlet pressure is set to 6.89 MPa. Liquid is injected at the bottom of the tube and heated over a length = 1 m. Table 6 sums up the inflow and wall boundary conditions specific to each of the three simulated experiments. Since the flow is axisymmetric, only a small angular section of the tube is discretized. The resulting mesh counts 16,800 hexahedral cells.
Cross section void fraction is measured over the tube length. In the heated lower section of the pipe, subcooled boiling occurs and steam is generated. The section above is adiabatic and vapor condensation occurs due to the mixing of the vapor generated near the heated wall in the lower section with the still subcooled liquid core. Figure 51 shows that new and standard model gives similar results and the results are in good agreement with experimental data.
Boiling Flow in a Vertical Pipe, with Refrigerant Fluid as Simulating Fluid. This experiment, carried out at the CEA (Commissariat Energie Atomique), is devoted to the study of upward boiling bubbly flows in a vertical pipe with a circular cross section [44]. The subcooled working fluid, Freon R12, is injected in the form of pure liquid at the bottom of the tube. The length of the pipe (5 m) can be divided into three successive parts: an adiabatic length (1 m), a heated length (3.5 m), and an adiabatic outlet section (0.5 m). The measurement section is located at the end of the heated length and the internal diameter of the tube is equal to 19.2 mm. We summarize the operating conditions for the selected tests in Table 7.
Because of the axisymmetrical nature of the experimental facility and the flow, the simulations are carried out in an angular sector of the section of the pipe. Meshes are generated as three-dimensional ones, with only one cell in the orthoradial direction, 20 cells in pipe cross section, and 200 cells along pipe axis. Let us note that mesh convergence has been achieved in previous works (Mimouni, 2010, 2011).
Figures 52 and 53 show that both models give similar results; the calculation of the diverse parameters of nucleate boiling (nucleation site density, diameter, and frequency) leads therefore to similar values with a different way.
Figure 54 shows that velocities are not modified by the new parameters.
In order to show that the wall temperature is better predicted with the new model, it is important to have experimental results on the wall temperature. Manon [44] has provided experimental data on wall temperature from DEBORA test. In Figure 55 heat flux, mass flow rate and pressure are fixed and only the inlet temperature is varying. Results show clearly that the wall temperature is overestimated with the standard model but the new model predicts accurately the wall temperature. This is of relevant interest for critical heat flux calculation.
Boiling Flow in a Vertical Annulus, with Refrigerant Fluid as Simulating Fluid. The Arizona State University (also known as ASU) experiment is described in [45, 46] for two-phase boiling flow measurements. The test section of the ASU experiment consists of a vertical annular channel with a heated inner wall and an insulated outer wall. The inner tube is made of 304 stainless steel (inlet diameter: 14.6 mm; outlet diameter: 15.9 mm) and the outer pipe of transparent pyrex glass (inlet diameter: 38.1 mm; outlet diameter: 47 mm) except for a 0.496 m long measurement section which is made of quartz (inlet diameter: 37.7 mm; outlet diameter: 41.7 mm). The inner tube is resistively heated, the upper 2.75 m of the 3.66 m long test section being the heated length. The lower 0.91 m serves as the hydrodynamic entrance length. The working fluid is a Freon R-113 (1,1,2-trichlorotrifluoroethane). We give the operating conditions for the selected test case in Table 8.
For the facility geometry, radial symmetry has been assumed, so that the numerical simulations could be performed on an angular sector of the pipe with symmetry boundary conditions at both sides. Three different grid sizes are used in this study (Table 9).
Standard and new model give similar results for void fraction and liquid temperature profiles in reasonable agreement with the experimental data (Figures 56-57).
This experiment is a bit special given that it is a low Reynolds number case. That is why the velocity profiles are not well predicted near the wall (Figure 58). In fact, the turbulent viscosity is too low in this configuration to have an accurate prediction. In NEPTUNE_CFD, the turbulent added term for -ε model may be underestimated; it would explain the overprediction. Moreover, a low Re turbulence model should be used instead of the classical SSG model used in the calculations. Further calculations will be performed with the -ε EBRSM model available in the code. Another consequence of the low Reynolds number can be observed in Figure 57 where the liquid temperature is underestimated with the finest grid. Indeed, the forced convective heat transfer coefficient given by Dittus-Boetler is not well adapted very near the wall where the velocity is the lowest.
Nevertheless, a reasonable agreement is reached on a case where it is not fitted for.
Boiling Flow in Rod-Bundle, with Refrigerant Fluid as Simulating Fluid. The objective of this test case is to verify that the new nucleate boiling model not only gives reasonable agreement in simple geometry where the original correlation is working well but remains still valid in more complex industrial applications where the geometries are not simple vertical tubes (e.g., PWR fuel assemblies).
POSEIDON is a rectangular channel with three heated pipes inside carried out at EDF R&D. The liquid is the Freon 12 (R12) and is going up around the tubes in order to simulate thermohydraulic conditions close to the real ones from pressurized water reactor (PWR). The scale of the experiment is twice bigger than that in reality; the diameter is 19 mm (8 mm in reality) with 4.25 mm between each of them. The heated length is 2 m. The pressure is 25 bar, the inlet temperature is 64°C, the inlet mass flux is 6 kg/s, and the heat flux is 238,810 W/m2 which is extremely close to the condition of use of PWR as the R12 is a simulating fluid. Therefore, the experiment is close from an industrial case. Calculations have been performed with three levels of mesh refinement (100,000, 500,000, and 1 million of cells).
Only few experimental results are available and measurement probes are represented in Figure 59. For three different meshing refinements, calculations have been carried out. Figure 60 shows that the liquid temperature agrees fairly with experimental data. As previous cases, the mesh convergence is acceptable.
2.6. Conclusion about Our Current Capabilities to Simulate Boiling Flows
In the framework of the nuclear industry, a CFD tool has been developed and advanced models dedicated to boiling flows have been implemented and validated against experimental data for ten years now including a wall law for boiling flows, wall heat transfer for nucleate boiling, and turbulence and a polydispersion model. A mechanistic modelling of CHF taking into account all the technical results and the experience gained was firstly proposed in [38] and generalized in the present paper. Calculations have been performed on 1500 experimental cases covering a large physics scope in terms of mass flux, pressure, quality, and channel diameter. Wall temperature excursion and onset of boiling crisis are reasonably well reproduced.
The results found are of central importance in a wide range of industrial applications and will be applied to industrial geometries in further studies.
But, the general tendency was to overestimate the wall temperature. As a consequence, a new nucleate boiling model based on the Liu-Winterton correlation has been validated for water and freon boiling test cases. The validation cases show that liquid temperature, gas velocity, and void fraction are similar for both nucleate boiling models. But the wall temperature is calculated more accurately in comparison with the standard model. Thus, the new nucleate boiling model improves the modelling of boiling flows in NEPTUNE_CFD.
Moreover, the range pressure validity is much wider for the new nucleate boiling which is of relevant interest in accidental scenarios. But this point is needed for supplementary validations.
The next step will be to combine the DNB model and the new nucleate boiling model proposed in this section. This work is in progress.
Moreover, in this part, we have only considered a uniform heat flux imposed at the wall. It can be considered as a first step towards industrial applications more related to nonuniform heat flux repartition. This important topic is currently assessed.
3. Modelling and Computation of Cavitation with a Two-Phase Flow Approach Downstream of an Orifice
Cavitation is a subject of interest in a number of fields in science and engineering [47]. The consequences are quite important in safety analysis of nuclear power and chemical and process plants. Cavitation may occur in numerous pipe components such as valves, orifices, elbows, and tees, in any place where low pressure and/or high velocity are/is reached. The effects of cavitation developments are first noise and vibrations, then erosion of the downstream pipe if the material is not resistant enough, and finally loss of performance of the component.
Most of numerical simulations of cavitation phenomena combine the Volume-of-Fluid (VOF) technique with a model predicting the growth and condensation of bubbles [48–51]. The bubble-liquid flow is treated as a single-fluid homogeneous model: balances of mass and momentum are solved for only one fluid. The flow can be supposed to be isotherm; therefore the energy balance is not solved. The flow is supposed to be incompressible.
The bubble growth is described by the Rayleigh-Plesset equation, applicable in the range of moderately low pressure. The liquid and vapor densities are most of the time given constant or are simply related to the pressure. The vapor pressure is given constant.
This model needs the nuclei concentration per unit volume of pure liquid, that is, the initial value of void fraction. From a physical point of view, vapor bubbles are supposed to be present in the flow independently of cavitation phenomena. Moreover, initial bubble diameter, characteristic times of growth, and condensation must be prescribed. In engineering applications, calculation parameters must be adapted to establish satisfying agreement with experiments, because of a large variety of cavitation types. Then these values are used to predict cavitation behaviour in other conditions.
This section deals with the numerical simulations of cavitation phenomena from a different point of view:(i)We solve, as already explained in previous sections, the same equations of the two-phase flow model (so-called 6-equation model): mass, momentum, and energy balance for both liquid and vapor.(ii)We solve the equations of energy to take into account compressibility effects.(iii)We solve the equations of momentum to take into account the drift velocity between bubbles and liquid.(iv)We solve an additional equation on the bubble diameter.(v)We solve compressible and unsteady flow transports equations.(vi)The thermodynamic properties of the real fluid are calculated everywhere in the flow and at each time step.(vii)Cavitation comes from preexisting cavitation germs (noncondensable gas and/or vapor bubbles and/or solid particles) in the flow.
The development of cavitation phenomena depends on many parameters, including the pressure, flow rate, depressurization rate, surface conditions, presence of dissolved gases or impurities, and local surface irregularities which may serve as nucleation and/or cavitation sites. In our model, microbubbles expand in regions where the local pressure is below the saturation pressure in the case of bubbly flow, with a tendency to agglomerate into slug bubbles.
This section is organized as follows. First we describe the set of equations solved in the numerical solver. Next, turbulent terms and interfacial transfer terms are detailed. The mass transfer modelling proposed in this section is detailed below. In the last part, the two-phase flow model is validated by simulating the EPOCA test.
3.1. Interfacial Transfer Terms
The cavitation inception may occur when cavitation nuclei in a liquid are exposed to a sufficiently low pressure, comparable with the vapor pressure of the fluid, to cause their unstable and explosive growth. The cavitation is believed to take place in low pressure regions, for example, produced by vortical structures where pressure decreases.
In the case of bubbly flow, the mass transfer for cavitating flow in most of engineering calculations is given by [52]The bubble diameter is calculated by solving an additional equation on the bubbles density in the flow. Initially, bubble diameter is supposed to be equal to few micrometers.
Moreover, cavitation phenomena are sensitive to vortical structures which are produced downstream of an orifice. The pressure may drop under the saturation pressure inside vortical structures. But, the overestimation of the turbulent viscosity with a EVM model is one of a major difficulty to capture vortical structures accurately. Another difficulty is the inability to control the grid refinement in vortical structures to calculate pressure inside vortices correctly. Therefore, the cavitation models account for the pressure drop in vortical structures by simply raising the phase-change threshold pressure value aswhere is saturation pressure corresponding to the local temperature.
In a similar way, Zhang et al. [52] take into account the turbulence-induced pressure fluctuations by where is the turbulent pressure. The modelling of is potentially of relevant interest for cavitation and needs to be carefully investigated in future calculations.
3.2. Computation of Cavitation Developments Downstream of an Orifice
The objective of this section is to compute the cavitation developments downstream of an orifice (Figure 61). The internal diameter is = 200 mm. They are mounted in the middle of a long straight pipe (42D long). The flow is assumed to be axisymmetric; therefore a two-dimensional axisymmetric meshing is used. We only compute 4.4 m length and the number of cells is about 10650. The meshing is refined near the orifice and the cells have an average length of 4.5 mm (Figure 62). The orifice has 4 mm thickness. The experimental tests have been performed using the EPOCA facility at EDF R&D in France [53]. The test conditions are controlled using two independent parameters: the flow rate and the pressure level.
Computations of cavitation developments downstream of a orifice have been performed with β = 0.4, where β = is the orifice opening coefficient.
The diameter of the vapor bubble is limited to 1 mm in the calculations. The geometry of the diaphragm is given in Table 10. The experimental conditions are given in Table 11.
Tests 15–18 contain nearly no vapor bubbles which agrees well with CFD calculations. Two-dimensional pressure fields clearly show low pressure zones (Figures 64–66). The cavitation level seems well predicted for tests 19–21 but the void fraction is overestimated just behind the diaphragm and near the wall. One of the major difficulties encountered during this study is the inability to evaluate the density of the cavitation germs because they mainly determine the order of magnitude of void fraction.
A core potential flow is observed just behind the diaphragm and centred on the axis. Cavitation phenomena are initiated outside this region. These results are in accordance with those quoted by [47].
The static pressure profiles (Figures 69–75) are in good agreement with the experimental data. A good prediction of the pressure is crucial to understanding the phenomena of noise and vibrations and then erosion of the downstream pipe.
3.3. Conclusion about Our Current Capabilities to Simulate Cavitation Phenomena
The simulations of the cavitation phenomena presented in this document consist in the resolution of balance of mass, momentum, and energy balance for both liquid and vapor phase. We consider turbulent, compressible, and unsteady flows. The thermodynamic properties are given by tabulated laws for real fluids. Moreover, cavitation inception occurs by preexisting cavitation germs in the flow. Microbubbles travel through the flow and grow in the regions where the local pressure is below the saturation. Void fraction condenses in the liquid region in a symmetrical way. The turbulence model has a predominant role because it describes the formation and behaviour of the vortex generated in the flow. The vortex breaks the slug bubbles into small bubbles. The small bubbles are captured by vortical structures, regions at low pressure, favourising vapor creation. Moreover, cavitation phenomena are sensitive to vortical structures where description is strongly related to the turbulence modelling. As a result, a second-order turbulence model based on a -ε approach for two-phase flow modelling should be used in further calculations.
The model predictions compared with the experimental data showed that satisfactory agreement could be obtained against experimental data of pressure. Only a satisfactory qualitative agreement is obtained for void fraction. The main learning is that validation against pressure profiles is not enough and measurement of void fraction profiles is required in order to understand the cavitation phenomena.
Many issues still require an important R&D effort as forces exerted on bubbles, bubble polydispersion, and the effect of noncondensable gases. Evaluating the density of the cavitation germs remains the main challenge because it determines the order of magnitude of void fraction.
4. Perspectives: Multifield Approach and Interface Locating Method for Two-Phase Flows in Nuclear Power Plant
Many situations in nuclear power plant are characterized by liquid-vapor interfaces such as boiling crisis and flows occurring in a steam generator. The industrial context is described in a few words in the following.
During a loss-of-coolant accident, a relevant safety issue in the context of PWR life extension study is the possibility of a cold shock, the so-called Pressurized Thermal Shock (PTS), when subcooled water flows from Emergency Core Cooling (ECC) system through the cold leg (a horizontal pipe) towards the reactor pressure vessel. The thermal-hydraulics problem is to evaluate the heating up of water between the ECC injection and the vessel. A two-phase stratified flow may take place in a part of the cold leg, as well as more complex turbulent two-phase flow configurations in the cold leg region close to the ECC injection and in the downcomer region below the junction with the cold leg. In the situations of interest, the cold leg is thought to be a liquid-vapor stratified flow; the ECC jet is thought to not break up and to plunge into the free surface; direct contact condensation takes place with a maximum in the ECC region of the cold leg; the free surfaces are thought to be mostly wavy or rough; the water and gas are turbulent.
In this context, the interface locating method is of relevant interest to deal with mass, momentum, and energy exchanges through the interface. In cells where the interface is not captured, the interfacial area, mass, momentum, and energy transfer terms correspond to those of bubbles, droplets, or other intermediate states.
Moreover, in two-phase CFD involving large interfaces, the liquid or gas volume fractions cover the whole range from 0 to 1 and the great majority of cells are of course out of large interfaces. There can be bubbles entrained below the jet or droplets ejected. These phenomena could require a multifield approach to be taken into account.
When a liquid is flowing onto a heated wall, the heat transferred by the wall to the liquid causes the liquid to partly evaporate, therefore giving a two-phase bubbly flow along the wall. This kind of situation, called nucleate boiling regime, is an efficient manner to evacuate the heat produced into the wall, for example, by Joule effect or by a nuclear reaction.
Two-phase flows are featuring many other industrial applications such as heat exchangers and chemical reactors. Based on the interface structures, several flow regimes can be identified and commonly separated into three main groups: separated flows, dispersed flows, and the last group that would contain the flows such as bubbly annular, churn turbulent, or slug flow. Although this classification has been experimentally confirmed since a few decades, the numerical simulation of complex two-phase flow regimes is still challenging and a universal model remains to be established.
There are several approaches for two-phase flow modelling, describing the interfaces with either a dispersed or a located point of view. Bubbly flows are often modelled with an Eulerian dispersed description, within the two-fluid model of Ishii [13]. The averaged momentum balance equation is in this case closed with a set of interfacial forces such as drag, lift, virtual mass, and turbulent dispersion. These forces rely on empirical or statistical correlations making assumptions on the bubbles mean diameter and shape, generally considered as spherical or slightly ellipsoidal. To accurately treat problems where these assumptions are not verified, an interface tracking model is preferred to follow the distortions of the bubbles.
On the other hand, large interface flows such as slugs or free surfaces are mostly simulated through located approaches such as front tracking [54], level-set [55, 56], or Volume of Fluid [57] with an Eulerian point of view or Lagrangian grid methods (Hyman et al., 1984) with a Lagrangian point of view. All these methods aim at calculating the local characteristics of the interface, such as the curvature and the normal vector, to model the interfacial transfer in the momentum equation.
New approaches are explored to simulate more accurately the transition regime between bubbly and separated flows. The concept of a four-field and two-fluid model has been presented and studied over the last decade [58, 59]. Multifield approach has been studied by treating the different phases as a single fluid with variable material properties and situations featuring interphase heat and mass transfer have been investigated by Lakehal et al. [60]. This is no longer the case in the paper, where mass, momentum, and energy balances are solved for each phase.
Each phase is decomposed into a continuous and a dispersed field, resulting in a four-field system of mass, momentum, and energy equations. This kind of approach requires the set of mass transfer terms between the continuous and the dispersed fields of the same physicochemical phase. A spatial cutting length is dividing a phase between unresolved structures that are modelled and the larger ones that are simulated. This concept allows the simulation of a wide range of two-phase flow regimes with both a good accuracy on the behaviour of the most distorted interfacial structures and less CPU consumption than the direct simulation of every two-phase scale. If the dispersed fields are commonly dealt with from an Eulerian point of view, several methods can be used to locate the interface between the liquid and gas continuous fields.
A hybrid multifield approach is based on this four-field concept. It consists in modelling the two-phase flow with an Eulerian approach, the gas phase being split in two fields. The small bubbles, assumed to be spherical, are modelled with a dispersed approach whereas the larger bubbles, considered as too distorted to be accurately described by correlations, are simulated through an interface locating method. As a first step towards this new approach, we have simplified the general concept by considering the liquid phase as continuous in our simulations. The main results obtained with NEPTUNE_CFD are presented in [61–64].
Finally, this multifield approach will be in future calculations an efficient tool to run all models (described in the paper and validated separately) together.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work has been performed in the framework of the NEPTUNE project, financially supported by CEA (Commissariat à l’Énergie Atomique et aux Energies Renouvelables), EDF, IRSN, and AREVA-NP, and in the framework of the MAGESTIC project, financially supported by EDF.