CFD Applications to Pressurized Thermal Shock-Related Phenomena

Abstract

In pressurized water reactor accident scenarios, the injection of water from the emergency core cooling system (ECCS) (ECC injection) might induce a pressurized thermal shock (PTS), affecting the reactor pressure vessel (RPV) integrity. Therefore, PTS is a vital research issue in reactor safety, and its analysis is essential for evaluating the integrity of RPVs, which determines the reactor life. The PTS analysis comprises a coupled analysis between thermal–hydraulic and structural analyses. The thermal–hydraulic approach is particularly crucial, and reliable computational fluid dynamic (CFD) simulations should play a vital role in the future because predicting the temperature gradient of the RPV wall requires data on the transient temperature distribution of the downcomer (DC). Since one-dimensional codes cannot predict the complex three-dimensional flow features during ECC injection, PTS is one reactor safety issue where CFD simulation can benefit from complement evaluations with thermal–hydraulic system analysis codes. This study reviewed from the viewpoint of the turbulence models most affecting PTS analysis based on papers published since 2010 on single- and two-phase flow CFD simulation for the experiment on PTS performed in the Rossendorf coolant mixing model (ROCOM), transient two-phase flow (TOPFLOW), upper plenum test facility (UPTF), and large-scale test facility (LSTF). The results revealed that in single-phase flow CFD simulation, where knowledge and experience are sufficient, various turbulence models have been considered, and many analyses using large eddy simulation (LES) have been reported. For two-phase flow analysis of air–water conditions, interface capturing/tracking methods were used in addition to two-fluid models. The standard and shear stress transport (SST) models were still in the validated phase, and various turbulence models have yet to be fully validated. In the two-phase flow analysis of steam–water conditions, many studies have used two-fluid models and Reynolds-averaged Navier-Stoke (RANS), and NEPTUNE_CFD, in particular, has been reported to show excellent prediction performance based on years of accumulated validation.

1. Introduction

When a hot reactor pressure vessel (RPV) wall cools rapidly, a heat load called a pressurized thermal shock (PTS) occurs, a critical nuclear safety issue [1]. In deterministic evaluation, PTS scenarios are analyzed for main steam line break (MSLB) accidents, large break loss-of-coolant accidents (LBLOCAs), and small break loss-of-coolant accidents (SBLOCAs) (JEAC4206-2007 [2]). In probabilistic evaluation, the stuck-open valve (SO) is analyzed in addition to the above three events. The scenarios contributing to penetration vary depending on the progress of material irradiation embrittlement [3], and the contribution of LBLOCA increases as irradiation embrittlement progresses. An example of a PTS scenario occurs when the emergency core cooling (ECC) system (ECCS) is activated during a loss-of-coolant accident (LOCA) in a pressurized water reactor (PWR), and cold water is injected into a cold leg (CL) filled with saturated water and steam (ECC injection). Then, a temperature stratification is formed due to the density difference of the fluids, and direct contact condensation (DCC) occurs at the interface. The stratified cooling water creates a large temperature gradient on the wall surface and generates strong thermal stress on the vessel’s inner surface because of stratified cooling water flowing from the CL to the downcomer (DC) and rapidly cooling the RPV wall surface. Therefore, if cracks (flaws) exist on the vessel’s inner surface and the thermal stress on the inner surface exceeds the limit, the crack can propagate and damage the vessel. In other words, since PTS degrades the integrity of RPVs and limits their operational life [4], the detailed thermal–hydraulic behavior during coolant mixing is essential for the integrity assessment of RPVs [5]. However, the SO is a phenomenon in which crack propagation occurs due to pressurization of the system as the relief valve closes after the wall has cooled (i.e., crack propagation occurs due to pressurization rather than thermal stress). Although the Japanese deterministic evaluation method does not require SO analysis, the US recommends considering SO risks.

Figure 1 shows the fundamental phenomena for thermal–hydraulics in PTS, as shown in the framework of the European Platform for Nuclear Reactor Simulations (NURESIM) project [5–8]. Some physical phenomena are still being elucidated, and the overall system simulation during the ECC injection process and the accurate reproduction of the RPV heat load present significant challenges [5, 9–11]. Therefore, it is necessary to quantitatively grasp the wall cooling from the CL to the DC. In the one-dimensional model, the complex mixing phenomenon of the DC is expressed by averaging, so the detailed simulation of the local phenomenon is ineligible [11, 12]. However, computational fluid dynamic (CFD) analysis can consider the details of the flow path geometry and predict multidimensional features of the mixing process between the hot water in the primary system and the ECC subcooled water or two-phase mixture in the CL and the DC. Consequently, the fluid temperature distribution in the DC has been concluded that the accuracy improved, and conservatives of the heat transfer coefficients on the inner RPV surface decreased [13, 14], and the CFD simulation for PTS is a central theme in many international projects. Numerous benchmark tests have been performed on separate or combined effect phenomena of PTS [5, 15–19]. Moreover, many experimental facilities have been built to provide high-resolution data for validating CFD codes and their physical models, and the focus has been on evaluating turbulence models against experimental databases [20]. This is because the effect of DCC occurring at the free surface of cooling water and stratified flow is the most important for predicting the temperature of the part subjected to thermal shock, and this process is strongly dependent on turbulence to promote mixing [5].

Figure 1 

Fundamental flow phenomenon during a PTS situation with a partially filled CL [7].

Therefore, this study is aimed at summarizing, from the viewpoint of the turbulence model, CFD simulations for single- and two-phase flows related to PTS that have been conducted. Reviewed reference papers were on validating the code and model by comparing the simulation results with the PTS experiment data. However, we do not mention the uncertainty quantification (UQ) analysis here because it has rarely been performed. Moreover, the simulations were categorized into turbulence models (RANS or LES) and fluid phases (single- or two-). Comparisons of results are only presented for common and representative results performed under the same experimental conditions, showing that particular models perform better under certain conditions. A review of two-phase flow CFD simulations up to 2007 has been summarized in detail by Lucas et al. [5]; therefore, this paper focuses on studies since then. For that purpose, the four representative large-scale experimental facilities for benchmarking CFD code for PTS are featured in Section 2. The CFD simulations on the single-phase flow to PTS were reviewed from the viewpoint of turbulence models in Section 3, and those on the two-phase flow were reviewed in Section 4. Finally, the review’s conclusions about the current situation of the CFD simulation studies to PTS are presented in Section 5.

2. Major PTS Experimental Projects

Demonstration sub- and full-scale experimental facilities for benchmarking CFD codes for PTS performed are as follows [5, 7]: for (low-pressure) single-phase flow, the Rossendorf coolant mixing model (ROCOM) test facility [21] and the FORTUM PTS mixing test facility [22], and for the two-phase flow (conditions for the single-phase flow can also be performed sometimes), HYBISCUS test facility [23], upper plenum test facility (UPTF) [24], large-scale test facility (LSTF) [25], condensation due to safety injection (COSI) test facility [26], transient two-phase flow (TOPFLOW) test facility [27], and experimental setup for the Organisation for Economic Co-operation and Development (OECD)-Texas A&M University (TAMU) CL mixture benchmark [28]. ROCOM, UPTF, LSTF, and TOPFLOW are large-scale experimental facilities described in detail below. Table 1 outlines large-scale experimental facilities.

Table 1 

Large-scale experimental facilities for PTS experimental projects.

2.1. ROCOM Test Facility (HZDR)

Helmholtz–Zentrum Dresden–Rossendorf (HZDR) conducted experiments in the ROCOM experimental facility [21, 29] (Figure 2(a)) to investigate the effect of primary loop inventory and ECC water density differences on DC mixing. The experimental facility comprises four loops and is a 1 : 5 linear scale of a German four-loop KONVOI-type reactor. When the loop volume is scaled to 1 : 125, the coolant transition time matches the actual reactor [30]. The core barrel with its lower support plate and core simulator, the perforated drum of the lower plenum, and the inlet/outlet nozzles are set inside the RPV. They are made of clear acrylic material and are intended for unheated experiments [31, 32]. DC velocity distribution is measured using a laser Doppler velocimeter (LDV), and the void fraction is measured at the reactor inlet, upper DC, lower DC, and core inlet using a wire-mesh sensors (WMS).

(a) ROCOM test facility [14]

(b) TOPFLOW-PTS test section [17]

(c) Primary system of UPTF [108]

(d) Schematic view of the LSTF

Figure 2 

Large-scale experimental facilities on PTS.

2.2. TOPFLOW Test Facility (HZDR)

The TOPFLOW test facility [27] of the HZDR is a multipurpose thermal–hydraulic facility investigating nuclear reactors’ steady- and transient-state two-phase flow behavior. The TOPFLOW-PTS test section (Figure 2(b)) was designed to verify the prediction accuracy of the two-phase flow CFD code on the two-phase flow behavior in the CL and DC during ECC injection. This experimental program reproduces the entire scenario of a PTS event under a small rupture LOCA [33]. The equipment can conduct two-phase steam–water experiments at a maximum temperature of 286°C and a maximum pressure of 7 MPa. It comprises the main reactor components in the pressure vessel: the flat DC, pump simulator (PS), ECC injection line, and CL [17, 34]. The experimental facility models an Électricité De France (EDF) CPY 900 MWe PWR at a scale of 1 : 2.5. The measurements used 196 thermocouples, a heat flux probe, local void probes with micro-thermocouples, an infrared camera with a frame rate of 10 f/s, a resolution of pixels, a measurement accuracy of approximately 1.5 K, a high-speed camera with a frame rate of 200 f/s and resolution of , and WMS of (the liquid velocity measurement based on the cross-correlation technique) [5, 33, 35]. This experimental device is constructed using diving chamber technology [33]. It is operated at pressure equilibrium between the test section and the internal pressure of the tank, and its design and structure are simpler, manufacturing costs are less expensive, and the thin-walled structures allow the use of specific measurement techniques.

2.3. UPTF (BMFT-KWU)

The UPTF [36] (Figure 2(c)) was erected and operated by the Siemens/Kraftwerk Union (KWU) by order of the Bundesministerium für Forschung und Technologie (BMFT). The UPTF aims to validate various ECCS concepts for nuclear power plants [37], and it is designed to study multidimensional two-phase flow effects, specifically in the upper plenum, upper core region, and DC [38]. The experimental setup is a 1 : 1 model of the four-loop 1300 MWe Siemens/KWU PWR in Grafenrheinfeld, Germany, with a primary system maximum pressure of 2.2 MPa. It fully mimics the PWR design, including the pressure vessel, upper plenum, DC, and core barrel. Therefore, there is no fidelity distortion of the two-phase flow simulation capability due to scaling. The measured parameters are the upper plenum and DC water level, mass flow rate, and thermocouple fluid temperature [4].

2.4. LSTF (JAEA)

Figure 2(d) is a diagram of the LSTF [25] of the Japan Atomic Energy Agency (JAEA). LSTF has been used in the OECD/Nuclear Energy Agency (NEA) Rig of Safety Assessment (ROSA) project, and the obtained result has been validated for computer code and simulation models. Steady-state single- and two-phase flow natural circulation experiments with the thermal stratification with subcooled water and the coolant mixing in the presence of superheated steam were performed, and three-dimensional temperature distributions of the CL and the DC during ECC injection were obtained [18]. The experimental equipment is a 1 : 48 scale of the 1,100 MW class Westinghouse-type PWR in actual volume, a 1 : 1 scale in actual height, 10 MW in electrical output, and 15.7 MPa in maximum operating pressure. The multidimensional thermal–hydraulic response was simulated using the reactor core comprising the 16 square fuel assemblies of approximately the same size as a PWR fuel assembly with 1008 (1/48 scale) sheathed electric heater rods and eight semicrescent fuel assemblies, an annular DC with a width of about 50 mm in a pressure vessel, and CLs with an inner diameter of 207 mm. Instrumentation has been set up with thermocouples, water level measurements, gamma-ray densitometers, and video probes to observe the flow visually.

3. Application of Single-Phase Flow CFD to PTS

Many PTS analyses for single-phase flows have been performed for experiments in the ROCOM test facility [21]. The single-phase flow CFD simulation conditions for experiments in the ROCOM test facility are shown in Table 2, the TOPFLOW test facility in Table 3, the UPTF in Table 4, and the LSTF in Table 5. The ROCOM test facility data are used besides PTS, and many CFD simulations on turbulent mixing have been performed [39]. However, this paper is limited to studies related to PTS. Moreover, Table 6 shows the characteristics of representative turbulence models that were used highly frequently in this paper. For single-phase flow, it has been reported that accurately reproducing the problem requires a turbulence model that accounts for low Reynolds effects, laminar-to-turbulent transition, and buoyancy effects [39].

Table 2 

CFD simulations of a single-phase flow on the ROCOM experiments.

Table 3 

CFD simulations of a single-phase flow on the TOPFLOW-PTS experiments.

Table 4 

CFD simulations of a single-phase flow on the UPTF experiments.

Table 5 

CFD simulations of a single-phase flow on the LSTF experiments.

Table 6 

Characteristic of representative turbulence models in this paper.

3.1. Application of Reynolds-Averaged Navier-Stokes (RANS) to PTS

Except for papers related to benchmarking analyses conducted under the International Atomic Energy Agency’s (IAEA’s) Cooperative Research Project (CRP) on applying CFD to nuclear power plant design, papers on CFD simulations using RANS on PTS for ROCOM data have been reported by Farkas et al. [40], Puragliesi [41], and Wei et al. [42] (Table 2). Both have been performed for experimental data conducted within the OECD/NEA Primärkreislauf 2 (PKL2) project [43].

Farkas et al. [40] investigated the effect of the differences in water density between the primary coolant loop and the ECCS on mixing in the DC. They targeted the experimental data in the ROCOM test 1.1 (12% density difference) in the OECD/NEA PKL2 project [43], applied 3D thermal–hydraulic CFD code of FLUENT [44] to create and optimize a numerical model (computational mesh), and verified the CFD simulation. A transient analysis was performed using the Reynolds stress model (RSM) [45–47] as a turbulence model. The calculated results qualitatively and quantitatively correlated with the experimental results and accurately predicted the experimental results, especially at the core inlet central region. However, a difference exists between the experimental and calculated results regarding the position of the minimum temperature at the wall side of the core inlet due to the various flow patterns in the DC and the corresponding temperature fields.

Puragliesi [41] performed numerical analyses for the ROCOM test 1.1 (12% density difference), test 2.1 (1.28% density difference) under quasi-steady flow conditions, and test 2.2 under transient conditions experiments in the OECD/NEA PKL2 project [43]. The density difference for test 2.2 equals test 1.1 but targeted the transient volumetric flow rate differing from test 1.1 volumetric flow conditions [41, 43, 48]. The applicability of the unsteady RANS (URANS) model has been validated using 3D thermal–hydraulic CFD code of STAR-CCM+ v10.04 [49] (the standard low-Reynolds (SLR) model formulated by Lien et al. [50] and the eddy diffusivity model) to evaluate the predictive capabilities of CFD codes for asymmetric flows in the natural circulation regime during an MSLB with the loss of off-site power accident scenario. Due to the increased turbulent diffusivity, the maximum/minimum value prediction of the mixing scalar (MS) was detrimentally affected, and the eddy diffusivity model showed problems for flows with high-density ratios and Froude numbers. They showed that scale-resolving simulation (SRS) should be used instead of URANS for estimating high-precision maximum and minimum local values and for specific scenarios with a decreasing Reynolds number with time. This is because of possible relaminarization in extensive regions of the RPV and a drastic change in the flow regime and structures.

Wei et al. [42] investigated the effects of buoyancy-influenced turbulent mixing in RPVs during MSLB accidents, targeting the experiments in tests 1.1, 2.1, and 2.2. They focused on the turbulent heat flux models’ effects intrinsically related to concentration/temperature, especially against buoyant turbulent mixing in large geometrical scales. They combined two turbulent models (the cubic nonlinear low-Reynolds model (LienCubicKE) and standard or buoyancy-modified shear stress transport (SST) model [51–54] and three turbulent heat flux models (simple gradient diffusion hypothesis (SGDH), general gradient diffusion hypothesis (GGDH), and algebraic flux models (AFMs)). The SGDH model assumes that the turbulent heat flux is proportional to the gradient of the Reynolds mean temperature field. It is suitable for forced convection. The GGDH model is ideal for shear-dominated flows because it introduces the dependence of the turbulent heat flux on the Reynolds stress tensor and has a certain level of anisotropy in the turbulent diffusivity. The AFM model includes a temperature dispersion term in the algebraic expression of the GGDH model’s turbulent heat flux and is suitable for strongly stratified natural convection [42, 55]. The 3D thermal–hydraulic CFD code used was OpenFOAM-5 [56]. Consequently, anisotropic turbulent heat flux models (GGDH and AFMs) showed improved prediction accuracy under strong buoyancy conditions and significantly improved the MS/temperature distribution prediction. However, in the case of weak buoyancy, they showed that the turbulent model plays a significant role, and the turbulent heat flux model has essentially negligible influence.

Among the CFD simulation for the ROCOM test in the OECD/NEA PKL2 project [42], we picked up the results of the quantitative comparisons on test 1.1 (12% density difference) performed by Puragliesi [41] and Wei et al. [42]. These results are about the spatially averaged MS at the DC and core inlet (CI) sensor sections (Figure 3(a)). From each calculation result, we selected the result closest to the experimental value (Puragliesi [41]: SLR model (turbulent Schmidt number ) and Wei et al. [42]: (standard) SST model and SGDH). Figure 3(b) shows the comparison of the results. The DC’s experimental result agreed with the SST model. The SLR model initially had a much better agreement in the CL, but the SST model was also within the error range. The mesh resolution of the calculation using the SST model is smaller, while enough accuracy is confirmed in the SST model.

(a) Sensor locations [41]

(b) Spatial-averaged MS at the position of DC (left) and CI (right) sensors [41, 42]

Figure 3 

Comparison with the experimental result on test 1.1 in ROCOM conducted in the OECD PKL2 project.

CFD simulation on UPTF data for single-phase flow was only analyzed using RANS. As shown in Table 4, Martin et al. [57], Beukelmann et al. [58], Höhne and Deendarlianto [59], and Li et al. [4] performed the analyses. Among them, Martin et al. [57] and Beukelmann et al. [58] have obtained boundary conditions using a thermal–hydraulic system analysis code, performed CFD simulation, and then fracture mechanic analysis.

Martin et al. [57] performed prediction accuracy verification and benchmarking among 3D thermal–hydraulic CFD codes of STAR-CD, Code_Saturne, and NEPTUNE_CFD (the latter two are coupled with the solid code SYRTHES) for single-phase flow PTS experiments based on UPTF. The turbulence model used a standard model. The pressure was 17 bar with an initial water and solid temperature of 190°C. ECC water was injected at a mass flow rate of 40 kg/s. The thermal–hydraulic system analysis code, the Code for Analysis of Thermalhydraulics during an Accident of Reactor and Safety Evaluation (CATHARE), was used to define the boundary conditions. Consequently, the calculation results of the three codes in the lower part of the CL agreed, and the DC fluid temperature correlated with all codes and experimental results. Therefore, they showed satisfactory calculation results in the three codes. Furthermore, a fluid–solid coupling was performed, showing that slight deviations between numerical computational and experimental results for the fluid did not affect the predicted solid temperature.

Beukelmann et al. [58] calculated pressure and temperature distributions using the thermal–hydraulic system analysis code Reactor Excursion and Leak Analysis Program (RELAP) and the spatial and temporal temperature distributions of the main circulation pipes and the RPV wall using the CFD code ANSYS CFX. The turbulence model used the SST model [60]. In the targeted UPTF-TRAM C1 series Run2a1, the RPV was filled with water, and the ECC water had a mass flow rate of 20 kg/s. The model used for these calculations was validated using the UPTF experimental results, and the model validity was confirmed by comparing the temperature distribution. Furthermore, a fracture mechanic analysis was performed using the obtained thermodynamic boundary conditions (such as fluid temperature at the inner RPV wall). Consequently, there is a large safety margin against brittle crack initiation for all investigated regions of the RPV and the most severe transients, and the RPV brittle fracture can be excluded.

Höhne and Deendarlianto [59] analyzed qualitative flow behavior using ANSYS CFX for UPTF test I for PTS. A standard model [61] was used as the turbulence model. In test I, the pressure was 1.8 MPa, and the fluid temperature, initial temperatures of the CL, RPV, and core barrel were 463 K (190°C). The ECC water at 300 K (27°C) was injected with a mass flow rate of 40 kg/s. Consequently, the stratification of the CL was predicted with satisfactory accuracy, and the calculated minimum temperature was within the experimental error. ECC water in the DC fell near vertically to the bottom of the RPV and had the coldest temperature. The flow behavior in the DC was well predicted, except for some nonphysical spurious circumferential oscillations.

Li et al. [4] considered including the strategy of meshing, time step, and turbulence model (standard [62], renormalization group (RNG) [63], standard [61], and SST [60] models) selections using Fluent 19.2 [64] to form best practice guidelines for CFD simulation of mixing phenomena in RPV. In the target UPTF test 11, the operating pressure was 1.8 MPa, the ECC injection temperature was 300 K, and the flow rates were 20, 40, and 70 kg/s. Consequently, the calculated results using the SST model correlated well with the experimental data, and the error was suppressed within ±5%. A noticeable local reverse flow occurred in the CLs and lower part of the DC. Furthermore, the reverse flow at the core entrance section improved the mixing effect. The reactor core cooling effect was better as the injection flow rate increased. However, as the flow rate increases, the fluctuation range of the pressure vessel wall temperature becomes larger, making it easier to cause thermal fatigue of the pressure vessel wall. Thus, it showed the need to pay more attention.

CFD simulation on LSTF for a single-phase flow only uses RANS. As shown in Table 5, Cai and Watanabe [65] performed numerical calculations on the preparatory test of test 1-1, and Farkas and Tóth [66] and Scheuerer and Weis [67] performed numerical calculations on test 1-1 [18]. As mentioned in Section 2, test 1-1 is an individual effect experiment (ST-NC-34) of LSTF conducted in the OECD/NEA ROSA project. The thermal stratification of the CL and DC during ECC injection under single-phase and two-phase natural circulation conditions is investigated (the two-phase flow analysis is described in Section 4.1), targeting quasi-steady-state problems at different mass inventories. In the experimental conditions, primary and secondary system pressures of a single-phase flow natural circulation were 10.28 MPa in the preliminary experiment and 15.5 and 6.7 MPa in test 1-1, respectively. The core power corresponded to 2% of the scaled nominal value, and the primary mass inventory was 100%.

Cai and Watanabe [65] used OpenFOAM-1.6 [68] to investigate the effects of the computational code, the turbulence model (standard [62, 69], RNG [70], and no turbulence models), the curvature of the CL between the pump exit downstream and the ECC injection nozzle, and the DC and pump exit on the temperature distribution. Reasonable results were obtained for the turbulence model with the standard model. The CL curvature indicated that the flow field was asymmetric when there was a curvature, and the high-speed region almost corresponds to the high-temperature region. In other words, they clarified that the curvature drives the mixing because the swirling flow is formed at the elbow. They also showed that the flow field and the mixing in the CL were more affected by the DC than at the pump outlet.

Farkas and Tóth [66] used FLUENT 6.3.26 [44] to investigate the effects, such as various turbulence models (standard , realizable models [71], and RSM) and spatial discretization schemes. The turbulence model affected the flow and temperature distribution in the CL, the standard model and RSM gave remarkably similar results, and the calculated results were significantly affected by buoyancy. In the realizable model, the cold plume was more intensively mixed, and the buoyancy effects were less dominant. However, CL thermal stratification was reproduced in all models and had little effect on the first/second-order discretization scheme. Moreover, the geometry of the CL and the DC connection part significantly affected the DC temperature distribution, indicating the importance of accurate geometrical modeling.

Scheuerer and Weis [67] used ANSYS CFX 12 [72] to investigate the effects of mesh resolution, discretization scheme, and turbulence model (SST model and Baseline- (BSL-) RSM [72]) on CL mixing. The SST model with low computational load was used for the turbulence model because no significant difference exists between the two models in the physical models’ investigation. Close to the ECC injection nozzle, turbulent mixing occurred at the shear layer between the CL water temperature and the ECC water, and the CL water temperature was underestimated. However, the temperature stratification in the CL and DC correlated well with the stratification data due to the buoyancy force and showed the superiority of the buoyancy effect.

Figure 4(b) shows the simulation result at the center of the cross-section of TE-2B4 and TE-3B4 (Figure 4(a)) on test 1-1 obtained by Farkas and Tóth [66] and Scheuerer and Weis [67]. Farkas and Tóth [66] used RSM, and Scheuerer and Weis [67] used the SST model. At the TE-2B4, the simulation result shows a similar tendency with no large difference between the models. Still, it does not agree with the experimental result, making predicting the turbulent mixing of ECC water injection a task. At the TE-3B4, the result of RSM is closer to the experimental result. We consider this because the turbulence mixing is promoted downstream, and the anisotropic turbulence becomes strong.

(a) Measurement plane locations [66]

(b) Temperature distribution in TE-2B4 (left) and TE-3B4 (right) [66, 67]

Figure 4 

Comparison with the experimental result on test 1-1 in LSTF conducted in the ROSA-V project.

3.2. Application of Large Eddy Simulation (LES) to PTS

Among studies on CFD simulation on PTS data acquired by ROCOM equipment, excluding papers related to IAEA CRP benchmark analysis, Loginov et al. [14, 31] and Feng et al. [73] published papers targeting LES (Table 2).

Loginov et al. [14, 31], under the FLOMIX-R project [19, 74], validated the turbulence model (LES subgrid model developed by Vreman [75]) for buoyancy-driven, transition, and momentum-driven experiments on the transient PTS performed in the ROCOM facility. The analysis target was the d05m00 test (5% density difference between a loop and ECC water and 0% flow rate in loop 1: buoyancy-driven) [31], one of 21 buoyancy-driven mixing experiments, and d05m05, 10, and 15 tests (5% density difference between a loop and ECC water and flow rates of 5%, 10%, and 15% in loop 1: buoyancy, transition, and momentum-driven) [14]. Buoyancy-driven mixing experiments investigated the effect of high-density ECC water mixing with low-density primary loop inventory water on the mixing and flow distribution in the DC. They varied the loop mass flow rate between 0% and 15% of the nominal flow rate, covering the range expected during the natural circulation mode, and they varied the density between 0% and 10% [16, 19, 74]. Under the same analysis conditions of 0% loop flow rate, the two simulations performed with slightly different turbulent fluctuations (turbulent disturbances) at the ECC inlet, and they were substantially different in DC mixing and quite sensitive to small turbulent disturbances in the ECC injection piping [14]. Similar phenomena were observed in the experiments, and the conventional comparison of a single experiment with the corresponding simulation could not directly quantify the LES accuracy for this flow. In the experiments with loop flow rates of 5% or more, they validated (1) the qualitative comparison of the single experimental and numerical realizations, (2) the application of ensemble averaging of experimental data, (3) spectral analysis, and (4) engineering approach to compare local temperature drops. The mixing phenomenon analysis correlated well with the experiment and correctly reproduced the flow regime in all regions considered. Regarding turbulent fluctuations, the power spectrum slightly overestimated 10% flow fluctuations but gave a conservative estimate of the temperature drop near the wall, showing a realistic and high reliability of the general CFD simulation.

Feng et al. [73] analyzed d10m05 (10% density difference between a loop and ECC water and 5% flow rate in loop 1: buoyancy-driven) of the previous buoyancy-driven mixing experiment using TrioCFD [76]. They compared it with the analysis conducted in 2005 using Trio_U for the same experimental conditions [77]. The turbulence model used the wall-adapting local eddy viscosity (WALE) model [78]. The dominant mixing phenomenon was well predicted, but overestimation of the experimental values was observed in the DC and perforated drum, and it was necessary to improve the computational mesh and address the high computational cost.

3.3. Application of RANS and LES to PTS

Under the CRP on applying the CFD cord to nuclear power plant design conducted in IAEA from 2012 to 2018, a CFD benchmark [16] for the d10m10 test was performed at the ROCOM facility (10% density difference between a loop and ECC water and 10% flow rate in loop 1: buoyancy-driven), one of the previously mentioned buoyancy-driven mixing experiments. This benchmark participated in eight institutions and used six main codes and six turbulence models. The final report [16] provided an overview of the projects and comparative results. Papers related to this benchmark include Höhne et al. [30], Höhne and Kliem [79], Čarija et al. [80], Chouhan et al. [32], and Ayad et al. [81], validating RANS or RANS and LES (Table 2).

Höhne et al. [16, 30] and Höhne and Kliem [79] performed simulations using ANSYS CFX and TrioCFD 1.7.3. The turbulence model used was RSM (BSL-RSM) for ANSYS CFX and the LES (WALE model) for TrioCFD (the computational mesh conditions used in RANS and LES differ). Furthermore, the code was compared under d00m15, which had not been examined until then (0% density difference between a loop and ECC water and 15% flow rate in loop 1: momentum-driven). Both codes and the turbulence model correlated qualitatively with the experimental data. The satisfactory reproduction of the transient dimensionless linear scalar quantity (MS) normalized according to the conductivity value obtained by the WMS (DC sensor) installed at the upper DC shows the capacity of CFD to analyze single-phase PTS events, even for the intermedia range between the momentum-driven and density-driven flows. The dominant mixing phenomena were predicted correctly. However, it was shown that the buoyancy forces in the mixed convection regime need to be calculated more accurately.

Čarija et al. [80] performed simulations using ANSYS Fluent [82] to investigate the effects of standard , standard , SST models, and RSM and the influence of various initial and boundary conditions on the solution. The SST model underestimated the concentration (temperature), while the standard model overestimated the concentration. The standard model was highly dependable (conservative) from a safety viewpoint. Other models overestimated the concentrations at some measurement locations but were invaluable for wall-bounded flows.

Chouhan et al. [16, 32] performed simulations using OpenFOAM-4.1 [83, 84] and the standard , SST models, and LES (one-equation eddy viscosity [85] and WALE turbulence models [78]). The computational mesh conditions used in RANS and LES differ. The LES predicted local and spatial mixing phenomena more accurately than the RANS-based model. The mean and maximum MS at all planes were better predicted using the WALE model than the one-equation eddy viscosity model. Compared with experimental data, the prediction accuracy of local mixing phenomena at all planes was satisfactory using a computational mesh with 12 M cells, the WALE model, the explicit second-order Euler method in time, and the central differencing in space.

Ayad et al. [81] performed simulations using ANSYS CFX 2022 R1 [86] and SST model. Furthermore, the Boussinesq approximation was used to consider density effects on the momentum equation. The mixing scalar distributions in the CL, DC, and core inlet were qualitatively and accurately predicted by SST compared to experimental data. However, the CFD results related to the turbulence phenomenon in the upper and lower DCs below the CL into which the ECC is injected are strongly underestimated. The physical phenomena that caused the quantitative deficit and the disturbance of the results were strongly related to the balance of frictional pressure that reduces the inertial effects of the direct contact between the inner wall and the mixture in the DC.

Among the IAEA CFD benchmark for the d10m10 test, Figure 5(b) shows the comparison of results for local MS at the CL center (at the position of 22.5° circumferential angle) (Figure 5(a)) of the upper DC sensor. If some turbulence models are tested in the paper, we selected the model with calculation results closest to the experimental result (Höhne et al. [16, 30] and Höhne and Kliem [79]: RSM with ANSYS CFX and LES (WALE model) with TrioCFD, Čarija et al. [80]: standard model with ANSYS Fluent, Chouhan et al. [16, 32]: LES (WALE model) with OpenFOAM, and Ayad et al. [81]: SST model with ANSYS CFX). RSM and LES, which consider anisotropy of the turbulence, yield a maximum value much closer to the experimental results than those obtained using the eddy viscosity model, demonstrating their superiority.

(a) Upper downcomer sensor locations [81]

(b) Local MS at the position of the center of the CL (circumferential angles of 22.5°) in the upper DC sensor [30, 32, 80, 81]

Figure 5 

Comparison with the experimental result on d10m10 in ROCOM conducted in the IAEA CRP.

Rodio and Bieder [87] and Bieder and Rodio [34] performed CFD simulation on TOPFLOW-PTS for a single-phase flow (Table 3). They performed the simulations on the temperature mixing of cold water injected from the ECC and saturated water in the CL of a 100% saturated water level of the CL diameter under high-pressure conditions of 5 MPa. In other words, in ROCOM, the flow is buoyancy-driven by the density difference of the isothermal-unheated liquid phase. However, in TOPFLOW-PTS, the flow is buoyancy-driven by the density difference depending on the temperature difference. The temperature difference was approximately 200 K, and the density difference was 200 kg/m³ (the density ratio was 0.2). Rodio and Bieder [87] used NEPTUNE_CFD 3.2.0 and TrioCFD. In the NEPTUNE_CFD simulation, the CL was always filled with the liquid phase by holding the water/vapor interface in the upper part of the DC approximately 5 mm above the CL. Therefore, since the condensation rate did not affect the thermal stratification progress in the CL, the flow is considered a single phase and analyzed using the two-fluid model code NEPTUNE_CFD. RSM (second-order Speziale–Sarkar–Gatski RSM ( SSG) [88]) was used for the transient analysis by NEPTUNE_CFD, and LES (WALE model) was used for the transient analysis by TrioCFD (the details of LES are described in [34]). This study compared the modeling methods for a compressible fluid being suitable for reproducing temperature and backflow, a dilatable fluid (low Mach number compressible fluid) that assumed that density is a function of temperature but not pressure, and an incompressible fluid by the Boussinesq approximation. The results indicated the effect of these models on buoyancy modeling in the CL. NEPTUNE_CFD used a compressible fluid model, and TrioCFD used incompressible (Boussinesq approximation) and dilatable fluid models. A comparison of the temperature distributions by the compressible and dilatable fluid models reproduced the experimental results well at all vertical positions. However, the incompressible fluid model showed more significant errors than the others. Moreover, although it was an analysis-only comparison, the axial velocity distributions correlated at all vertical positions, regardless of the modeling method. LES and URANS showed good overall agreement.

4. Application of the Two-Phase Flow CFD to PTS

Many PTS analyses for the two-phase flow have been performed targeting TOPFLOW-PTS experiments. So far, analyses have been performed addressing air–water, steam–water, and both conditions. Tables 7 and 8 show them classified into air–water and steam–water conditions. In two-phase flow, three types of turbulence have been identified near the interface for two-phase flows: turbulence diffusing from the wall boundary, turbulence generated by interfacial friction, and turbulence caused by interfacial waves. Thus, near the interface, it is necessary to consider the anisotropy of the turbulence [5]. The research on two-phase turbulent flow is still in progress, and one or more standard , , and SST models, RSM, and LES are used in the reviewed papers. It is assumed that RSM and LES were applied to consider the anisotropy of the turbulence near the interface, but eddy viscosity models are also often used. Most experiments were conducted under the European Union Nuclear Reactor Integrated Simulation Project (NURISP) framework [89, 90]. Here, the reference conditions for steady-state air–water and steam–water experiments [15, 91] performed on the TOPFLOW test facility were defined to improve the CFD modeling of two-phase flow PTS scenarios. However, the steam–water experiment could not be performed under this condition; therefore, no comparison was made with the experiment and only a pretest [15]. In the air–water experiment under this reference condition, the working pressure was 2.25 MPa, the CL water level was 50% of the CL diameter, and the ratio of the PS injection inlet mass flow rate to the ECC injection mass flow rate was 1 : 1.7. The temperature of the ECC injection was below the saturation pressure of the working pressure. The DC outlet mass flow rate is . The steady conditions were obtained with a constant water level in the CL. The initial water and air temperature in the CL was defined as a mixed temperature. In the steam–water experiment, the working pressure was 5 MPa, the CL water level was 50%, and the ratio of the mass flow rate of the PS injection to was 1 : 1.7. The PS outlet flow rate was to avoid condensation in the PS. However, the injection temperature of the PS fell below the saturation temperature, causing more steam to condense in the PS. Moreover, the water level in the CL must be constant to maintain a steady state, and the outlet mass flow rate of DC is (: total condensation amount). The steam temperature was taken as the saturation temperature, and the water and steam initial temperature in the CL was the saturation temperature. Detailed values, such as temperature, are given below [91]. The papers covering this reference condition on the TOPFLOW-PTS analyses are all performed in air–water (Table 7), and only the paper published by Apanasevich et al. [15, 91, 92] and Martin et al. [93] are performed in steam–water (Table 8) conditions. Furthermore, Table 9 shows UPTF-TRAM C2 Run6b in UPTF, and Table 10 shows the two-phase flow CFD simulation conditions for the preparatory test for the OECD/NEA ROSA project test 1-1 (also shown in Section 3.1) in LSTF.

Table 7 

CFD simulations of a two-phase (air–water) flow on the TOPFLOW-PTS experiments.

Table 8 

CFD simulations of a two-phase (steam–water) flow on the TOPFLOW-PTS experiments.

Table 9 

CFD simulations of a two-phase (steam–water) flow on the UPTF experiments.

Table 10 

CFD simulations of a two-phase (steam–water) flow on the LSTF experiments.

4.1. Application of RANS to PTS

Deendarlianto et al. [94], Deshpande et al. [95], and Kim et al. [96] performed the simulations targeting air–water conditions of the RANS-only analyses (Table 7). These referenced the experiment on the NURISP framework described above.

Deendarlianto et al. [94] demonstrated the feasibility of the Euler–Euler approach using the algebraic interfacial area density (AIAD) model [97–99], where the momentum exchange coefficient depends on the local morphology because the frictional resistance models of the air–liquid interface were required for the continuous phase with the dispersed and separated continuous phases, respectively. They also used the ANSYS CFX, implementing a new drag coefficient model. The two-phase flow approach used the two-fluid model with the AIAD model and the two-fluid model without interfacial momentum transfer (the standard two-fluid model). The fluid properties were based on the steam table of the International Association for the Properties of Water and Steam (IAPWS), and the turbulence model used the SST model. The temperature distributions in the CL and the DC obtained by the calculation correlated with the experimental values within the measurement error range. No temperature stratification was observed in the CL and DC, except for immediate upstream from the ECC injection, confirming that it is equal to the perfect mixing temperature.

Deshpande et al. [95] and Kim et al. [96] used Fluent 14.5 [100], the volume of fluid (VOF) method [100] as the two-phase flow approach, and the SST model as the turbulence model.

Deshpande et al. [95] employed three interface sharpening functions (compressive scheme, level set (LS) method, and free-surface model) [101] for the VOF method to improve the prediction performance for a third-order discretization scheme quadratic upstream interpolation for convective kinetics (QUICK). They compared the effect of the turbulence model using a standard model and the experimental result. Adding the interface sharpening functions improved the prediction accuracy of the fluid temperature distribution, and the free-surface model showed the best prediction accuracy. In the turbulence model comparison, both models accurately predicted the temperature distribution. The SST model showed better prediction performance than the standard model close to the injection.

Kim et al. [96] performed calculations using the VOF method to investigate the effect of the computational mesh resolution, evaluate the effect of computational mesh refinement in the PS and ECC injection region, and perform a sensitivity analysis of the turbulence model. The turbulence model used was the SST , and the standard model was used for comparison. The SST model significantly improved the temperature prediction near the wall. However, fine-tuning the bridging parameters (generally known as a blending function) is required at the gas–liquid interface. Therefore, the standard model remains dominant in the interfacial region.

Among the RANS-only analyses, Apanasevich et al. [91] and Martin et al. [93] performed pretest simulations in air–water (Table 7) and steam–water conditions (Table 8). The target experiments are the reference conditions shown above.

Apanasevich et al. [91] investigated the effects of the thermal mixing process using the homogeneous and two-fluid models for the two-phase flow approach in ANSYS CFX 12.0 [72]. Since this is a pretest simulation, it did not include a comparison with experimental data. Moreover, the heat transfer effect between structures and fluid was not considered. The turbulent flow model adopted the SST model. Regarding heat transfer, in the air–water condition, the homogeneous heat transfer model was used without considering the heat transfer in the air–liquid interface. The steam was assumed to be isothermal in the steam–water condition, and only one energy equation was applied to the water. The DCC model based on the surface renewal theory (model) proposed by Hughes and Duffey [102] was adopted for the latter condensation heat transfer between gas and liquid. Furthermore, to define the heat and mass transfer between the two phases, they used the two-resistance model built into CFX, which solves the heat transfer between each phase and the interface separately. Water, air, and steam properties were assumed to be constant. Consequently, under the air–water condition, the temperature distribution significantly differed in the region close to the ECC injection because the homogeneous and two-fluid models had different predictions of the jet split in the CL after the ECC injection. Under steam–water conditions, both models showed thermal stratification in the CL and the DC. The calculated total condensation rate using the two-fluid model was larger than the homogeneous model, and the difference was due to the cold-water plume formation and propagation prediction in the DC.

Martin et al. [93] performed pretest simulations using NEPTUNE_CFD to improve knowledge and understanding of crucial physical phenomena at the mock-up scale. They performed the simulations with a thermal coupling between fluid and solid domains. A two-fluid model was applied for the two-phase flow approach, and a standard model [69] was used for the turbulence model. The free-surface modeling (interfacial momentum transfer) used the large interface model (LIM) [103]. For heat and mass transport, there was a description that interphase heat transfer is modeled under steam–water conditions, but no specific model was specified. The air–water characteristics were constant, and the steam–water characteristics used the CATHARE steam table. Consequently, gas–liquid flow separation occurred in the ECC injection line when the flow rate fell below a specific threshold, significantly influencing the mixing phenomena in the impinging jet and downstream stratified CL area. Increasing the ECC flow rate confirmed the presence of a much higher hysteresis phenomenon for flow separation in the CL in the air–water condition than in the steam–water condition. Thermal stratification is formed upstream of the CL in the air–water condition. However, the air–water condition could not maintain the steady-state downstream and mixed perfectly, unlike the steam–water condition that can maintain steady-state stratification due to the DCC of the steam.

Apanasevich et al. [35, 92], Coste and Lecomte [104], Coste and Mérigoux [105] (overlap content with Mérigoux et al. [106]), Mérigoux et al. [17, 107] (Table 8), Cremer et al. [108] (Table 9), and Cai and Watanabe [65] (Table 10) performed analyses targeting steam–water conditions.

Apanasevich et al. [92] used ANSYS CFX 12.0 [72] to analyze the reference conditions of the NURISP project. They used the two-fluid model for the two-phase flow approach with and without (the standard two-fluid model) the AIAD model for interfacial momentum transport, the DCC model based on the surface renewal theory proposed by Hughes and Duffey [102] and the two-resistance model for thermal mass transport, and the IAPWS for fluid properties. The turbulence model was the SST model. The total condensation rate calculated using the AIAD model was higher than that using the standard two-fluid model. Furthermore, temperature stratification was observed at the entrance to the DC in the CL and the DC itself, indicating a large buoyancy effect. However, the cold-water plume formation significantly differed between the AIAD and standard two-fluid models. Accurate prediction of cold-water plume formation and the lack of formation was essential for structural analysis and general safety assessment and required further validation.

Coste and Lecomte [104] used TOPFLOW-PTS steady-state steam–water test (SSSW) 3-2 experimental data, validated in NEPTUNE_CFD 1.3 using the model (version 1.0.8) where LIM (excluding turbulence) is implemented as free-surface modeling with satisfactory results in previous validations [109]. They compared the results using an alternative condensation model (version 2.0-beta) in the steam–water interface. Details of the SSSW 3-2 test were not disclosed, but the water level was 50% of the CL diameter, saturated steam was supplied through a top pipe of the DC, and steam surplus was discharged to the condenser through an opening at the top of the DC. The two-phase flow approach used the two-fluid model, the free-surface modeling used the LIM, and the thermal mass transport modeling used the DCC model based on the surface renewal theory proposed by Coste [110], Lakehal et al. [111], and Magnaudet and Calmet [112]. The turbulence model used the standard model. The result with version 1.0.8 overestimated the heat and mass transfers and did not correlate with the experimental result, but the result with version 2.0-beta correlated well with the experimental data.

Subsequently, validation studies on the TOPFLOW-PTS database for steam–water data conducted in the NURISP project framework were taken over by the European Union Nuclear Reactor Safety Simulation Platform (NURESAFE) project framework [113, 114]. Under this project, Apanasevich et al. [35], Coste and Mérigoux [105], and Mérigoux et al. [17] used TOPFLOW-PTS steady-state steam–water experimental results to verify the reliability of the two-phase turbulent stratified-flow simulation with DCC using the two-fluid model. Although the paper by Apanasevich et al. [35] does not show the details of the experiment, it is presumed to be the targeting SSSW 3-17 because it corresponds to the CFX analysis result by Mérigoux et al. [17]. SSSW 3-17 has no water injection and discharge on the PS side, unlike the NURISP reference conditions. Therefore, the liquid phase flow rate at the DC outlet was controlled to maintain the CL water level at 50% of the CL diameter (steady-state). The constant reference pressure in the vapor phase was 5 MPa, and the pressure field was initialized with hydrostatic pressure. The ECC injection mass flow rate was , and the ECC water temperature was . The outlet mass flow rate on the DC is defined as ( is the integral condensation rate). The integral heat loss was also measured in the experiments but was not considered in the calculations because the effect on the condensation rate and temperature distribution due to heat loss was small [104]. The steam mass flow rate is , and the steam temperature equaled the saturation temperature (537.15 K), defined as the steam inlet and initial temperatures [17].

Apanasevich et al. [35] used ANSYS CFX 14.5 [115] for the two-fluid model with the two-phase flow approach, the AIAD model and the constant drag coefficient for interfacial momentum transport modeling, and two DCC models by Hughes and Duffey (HD model) [102] and Magnaudet and Calmet (MC model) [112] based on Higbie’s surface renewal theory [116] for interfacial heat and mass transport modeling. The turbulence model is the standard model. The simulation evaluated the modeling of interfacial momentum and heat transfer (the modeling water heat transfer coefficients determining interfacial heat and mass transfer due to DCC and the modeling drag coefficients representing interfacial momentum exchange between phases). The simulations using constant and variable drag coefficients show significant differences in results, confirming that the drag coefficient is a variable quantity that depends on the local fluid parameters. Moreover, the HD model overestimated the total condensation rate by a factor of two. The MC model also overestimated the CL and DC temperatures, but the total condensation rate could be predicted with experimental accuracy. The temperature distribution correlated well with the experiment. Therefore, it was concluded that the AIAD and MC models could be applied to predict temperature distribution and condensation amount during two-phase flow PTS scenarios.

Coste and Mérigoux [105] used NEPTUNE_CFD to compare with SSSW 3-16, 3-18, and 3-19 in addition to the SSSW 3-17 experiment shown above (the liquid flow ratios differed: SSSW 3-16, 3-17, 3-18, and 3-19 were , 1, 0.7057, and 0.4085, respectively). They investigated the ECC mass flow effect, compared the standard model and the new version with the SSG model suited for highly anisotropic turbulence, and mentioned the improved effect. The two-phase flow approach was the two-fluid model, and the free-surface modeling was LIM. The turbulence and DCC models were distinguished by the version used: version 2.0.0 was the standard model and the MC model [112], and version 2.2.0 was the RSM ( SSG) and the DCC model by Lakehal et al. [111]. The results were satisfactory, except for SSSW 3-18, where the ECC jet structure turned into stratified liquid vapor in the experiment [101] because in the analysis of SSSW 3-18, a stratified flow, such as a churn or slug flow, occurs, and the condensation rate fluctuates; therefore, it was challenging to predict the ECC flow progress. Also, version 2.0.0 overestimated the condensation rate and temperature, but version 2.2.0 significantly reduced the overestimation.

Mérigoux et al. [17] summarized the CFD benchmarks of the NURESAFE project framework, comparing calculations with ANSYS CFX 16.0 [117], NEPTUNE_CFD 3.0, and TransAT 5.1.2_RC2. Figure 6(b) compares the temperature distributions obtained by each code to the experimental results at the locations shown in Figure 6(a). CFX and NEPTUNE_CFD targeted the same experiments as previous papers [35, 105], which validated the turbulence model and computational conditions. However, CFX used a DCC model by Lakehal et al. [111], and NEPTUNE_CFD used a different version. Considering the SSG model in NEPTUNE_CFD 3.0 was performed by Mérigoux et al. [106] described later, it seems that the results [17] obtained here were reflected; however, the relationship with the NURESAFE project is not described. TransAT used the LS method of the one-fluid model for the two-phase flow approach and the model proposed by Coste [110] for the phase change model. Here, the mass evaporation rate based on the surface renewal theory is shown in terms of turbulent Reynolds and the molecular Prandtl numbers. The standard model was used as the turbulence model, which had better numerical stability than the very LES (VLES) (hybrid model of RANS and LES). In the CL, the NEPTUNE_CFD and CFX results similarly predicted the water temperature near the wall, and the NEPTUNE_CFD result correlated well with the experimental near the interface. However, CFX overpredicted the turbulent kinetic energy and eddy dissipation, underpredicting the water temperatures in the interface region and requiring improvements in the turbulence model. TransAT significantly underestimated the upstream and downstream water temperatures, but the DC inlet correlated the best with the experimental results. CFX and NEPTUNE_CFD correctly predicted the DC temperatures. TransAT failed to model the transport of cold ECC water correctly and, thus, significantly overestimated the water temperature. Reasons for this result could be a wrong prediction of the ECC jet behavior and the influence of buoyancy effects that are not considered using constant fluid properties (for numerical stability reasons).

(a) Thermocouple lance locations in the CL (left) and the DC (right) [15]

(b) Normalized temperature distribution in the CL (upper) and the DC (lower) [17]

Figure 6 

Comparison with the experimental result on the SSSW 3-17 case in TOPFLOW-PTS conducted in the NURESAFE project.

Mérigoux et al. [106] focused on the integral validations of the NEPTUNE_CFD using the RSM ( SSG model) for turbulence prediction and LIM to present the validation and application of NEPTUNE_CFD 3.0 under PTS phenomena. First, they confirmed the lattice dependence under SSSW 3-17 conditions. It was confirmed that it is within the uncertainty range of the experimental condensation rate, having been estimated with a precision of ±25%, regardless of the mesh refinement. Next, the conjugate heat transfer between liquid and solid was considered by coupling with the thermal–solid code SYRTHES developed by EDF, and the code was verified in the transient scenario transient steam–water (TSW) 3-5, including conjugate heat transfer. The water level in TSW 3-5 decreased until the CL was empty; the ECC flow temperature slowly reduced. Later, the water level rose again and stabilized when the CL was filled with water. Moreover, throughout the scenario, a rectangular steel plate instrumented with thermocouples contacted the water flowing in the DC. Due to the simulation, there was uncertainty about the boundary conditions to be applied to the solid part that was not in contact with the water to model heat conduction in the solid. From the experimental temperature profiles near the boundaries, they could deduce that some heat exchange occurred between the steel plate and its surroundings, but they could not quantify it. If they did not take into account these thermocouples closer to the boundaries, the maximum error of the steel temperature prediction decreased by up to 4%.

Cremer et al. [108] performed analysis for the UPTF-TRAM C2 Run6b experiments for the development and validation of applying the free-surface condensation model [102] based on surface update theory to ANSYS CFX (Table 9). The turbulence model used the standard model. In this experiment, the ECC water mass flow rate was 50 kg/s, and the DC water level was 2.7 m (measured from the bottom of the RPV test vessel), which allowed the formation of water streaks below the CL nozzle during ECC injection, where condensation occurred on the subcooled water surface. This result validated the free-surface condensation model and the drag coefficient for obtaining the maximum mass flow rate for cold-water attachment on the RPV wall. The drag coefficient obtained here was applied to the PTS transient analysis of the Gösgen–Däniken nuclear power plant (Kernkraftwerk Gösgen (KKG)) and compared with the existing KWU-MIX (analytical fluid-mixing codes validated with experiments) results.

Table 10 shows that Cai and Watanabe [65] used OpenFOAM-1.6 [68] to analyze the experimental conditions under the preparatory test for two-phase flow conditions of test 1-1 conducted at LSTF during ROSA-V program and verified the effects of the condensation model. The primary and secondary system pressure in this experiment was approximately 6.7 MPa. The injection mass flow rate from the ECCS was 0.3 and 1.0 kg/s for CLs A and B, respectively, and an additional water injection at a mass flow rate of 1.8 kg/s was performed for CL B under two-phase flow conditions. The two-phase flow model was the VOF method, the turbulence model is the standard model, and the condensation model implemented the simplified model comprising an interfacial area and heat transfer coefficient proposed by Watanabe and Nakamura [118]. The temperature distribution correlates well with the experimental results, confirming the effect of the condensation model.

4.2. Application of RANS and LES to PTS

Among the RANS and LES analyses, Ničeno et al. [119] and Apanasevich et al. [15] performed the analysis targeting the air–water reference conditions of the NURISP project (Table 7), and Apanasevich et al. [15] also performed the analyses targeting the steam–water reference conditions (Table 8).

Ničeno et al. [119] performed a transient analysis for air–water reference conditions using ANSYS FLUENT 12.0 [120] to investigate the effects of different turbulence models. The two-phase flow approach used was the VOF method, and the turbulence models were the standard , SST model, and LES (Smagorinsky and dynamic subgrid-scale (SGS) models [120]). Compared with the experimental results, the turbulence model has a considerable effect, and the temperature distribution using the standard model showed the best result. The simulation (steady-state analysis) using RANS by Deshpande et al. [95] shown in Section 4.1 also uses Fluent, the VOF method, the standard , and SST models in the same mesh resolution. However, the steady-state analysis showed better results than the transient analysis, and the SST model showed results closer to the experimental results than the standard model (Figure 7(c)). The physical time of transient analysis is 300 s, and sufficient convergence of the solution has been obtained [95]. The cause of this paradoxical result is unknown because no detailed information and consideration are presented.

(a) Thermocouple lance locations in the CL (left) and the DC (right) [15]

(b) Normalized temperature distribution in the CL (left) and the DC (right) [15]

(c) Normalized temperature distribution at the line LA1 using the VOF method with ANSYS FLUENT (Fluent)

Figure 7 

Comparison with the experimental result on posttest for the steady-state air–water case in TOPFLOW-PTS conducted in the NURISP project.

Apanasevich et al. [15] performed benchmark analysis using ANSYS CFX 12.0 [72], NEPTUNE_CFD 1.0.8, and ANSYS FLUENT 12.0 [120] to validate the CFD model for the two-phase flow PTS scenario. The papers listed above [92–94, 119] are highly relevant to this paper. Figure 7(b) compares the temperature distributions obtained by each code to the experimental results at the locations shown in Figure 7(a). Two-phase flow approaches and turbulence models used the two-fluid models and SST model in CFX, the two-fluid model and standard model in NEPTUNE_CFD, and the VOF method and LES (Smagorinsky SGS model) in FLUENT. CFX used AIAD and the constant drag coefficient as the interfacial mass transport, and NEPTUNE_CFD used LIM. The heat and mass transport in CFX used the HD model [102] and considered the heat transfer between the two phases. NEPTUNE_CFD used the DCC model [103] proposed by Lakehal et al. [121] based on the surface divergence model proposed by Banerjee et al. [122] and the DCC model [103] proposed by Coste and Lecomte [104] based on the surface renewal theory using the wall law to model heat and mass transfer by condensation. FLUENT used the HD model [102]. For the fluid properties, the air was regarded as a perfect ideal gas under the air–water condition of CFX, and IAPWS was applied to the liquid phase and the steam–water condition. They assumed that the air properties for the NEPTUNE_CFD air–water condition were constant (pressure 2.25 MPa, isothermal) and did not model heat transfer between phases. They applied the CATHARE steam table to define the liquid phase properties and the CATHARE steam table for the steam–water condition. The FLUENT air–water conditions were constant in each phase, the steam–water conditions did not use an equation of state assuming noncondensable gas, and the steam conditions used a user-defined function. For air–water conditions, comparisons were made with the experimental results, but for steam–water conditions, no comparisons were made because an experiment on reference conditions could not be performed. Therefore, under air–water conditions, FLUENT has the lowest temperature near the interface and the wall, and the temperature distribution near the interface correlates well with the experimental results. Figure 7(c) shows the verification result using Fluent [95, 96, 119], not included in Figure 7(b), among the analysis results from the same project shown in Section 4.1. Here, we organized the results using the VOF method, standard , and SST models. The results for the SST model with 5.7 M cells are the closest compared to those shown in Figure 7(b). The standard model with the same mesh resolution also agrees better than the SST model near the interface and is closer than the results for other conditions shown in Figure 7(b). It is seen that the influence of the mesh resolution is significant from this result.

Under steam–water conditions, RANS was more suitable to model two-phase flows with smooth and large free surfaces, indicating that the computational mesh size was unsuitable (too coarse) for LES. Furthermore, NEPTUNE_CFD had the highest condensation rate, whereas CFX and FLUENT had vastly different condensation rates despite using the same condensation model. They showed this result due to the impact of turbulent models on DCC modeling because condensation depends on the turbulence in the liquid phase: turbulent eddies transport the heat from the interface. Furthermore, only NEPTUNE_CFD and CFX predicted temperature stratification at the DC entrance, NEPTUNE_CFD predicted a uniform temperature very close to the saturation temperature at the DC, and CFX predicted a much lower nonuniform temperature than the other codes.

5. Conclusions

CFD simulations of PTS performed since 2010 were reviewed from the viewpoint of turbulence models, especially CFD simulations for the experiments performed in large-scale experiment facilities ROCOM, TOPFLOW-PTS, UPTF, and LSTF. The primary computing platforms were ANSYS CFX, ANSYS FLUENT (Fluent: from version 14.0) (formerly FLUENT), STAR-CCM+, NEPTUNE_CFD (focusing on the single-phase flow part under two-phase flow conditions), TrioCFD, TransAT, OpenFOAM, and Code_Saturne for single-phase flows and ANSYS CFX, ANSYS FLUENT (Fluent), NEPTUNE_CFD, and OpenFOAM for two-phase flows. Table 11 classifies papers’ target phenomena discussed in Sections 3 and 4 into single phenomena in the two-phase flow PTS situation listed by Lucas et al. [5] (including even research on single-phase flow). Phenomena not targeted in the paper but debatable based on the results were included in the phenomenon category at our discretion. The simulation’s main objects were turbulent mixing for single-phase flows and thermal stratification, interfacial momentum transport, and heat and mass transport for two-phase flows. In addition, it was also highly focused on heat transfer to turbulence production below the jet and at the wall in CL and heat transfer to the wall in DC, as shown in Table 11.

Table 11 

Reference classification based on single phenomena in the two-phase PTS situation [5].

In single-phase flow CFD simulation, where knowledge and experience are sufficient and reliable codes have been developed, turbulence models on LES, which are available to consider anisotropic turbulence, are increasingly being studied in addition to RANS. In the RANS, although the appropriate turbulence model depends on the flow pattern at the measurement position, the SST model and RSM were clarified to show better results than the other models overall from comparing the results. In the LES, many studies using the WALE model were especially reported, as the WALE model showed the best results among LES (e.g., [32]). Moreover, even with RANS, many analyses use computational meshes with a scale exceeding 1 M cells (elements or nodes). Table 12 shows the calculation time in the simulations targeted in this paper as a reference. It is expected that large-scale simulations will increase in the future as computer performance improves. At the same time, it is necessary to consider its use in industry, and it is essential to make an effort to reduce calculation costs while improving calculation accuracy.

Table 12 

Calculation characteristics.

In the two-phase flow analysis of air–water conditions, although it is necessary to consider the anisotropy of the turbulence near the interface, RANS based on the eddy viscosity model is still often used. In the two-phase flow analysis method, the interface capturing/tracking method (VOF and LS methods) analyses have been reported in addition to the two-fluid model. A comparison of the results shows that the result using the SST model and VOF method showed overall good agreement results in enough mesh resolution. When using the interface capturing/tracking method, the influence of the computational mesh resolution appeared significantly [96]. Analysis using the VOF method and LES did not show results superior to RANS’s because the computational mesh resolution is insufficient [15]. There is room for further investigation. In addition, in this study, we focused only on the fluid analysis by CFD, but in the simulation for UPTF, the thermal–hydraulic system analysis codes such as RELAP and CATHARE to determine boundary conditions and the coupling analysis with the structural mechanics (stress and fracture) analysis code were also performed [57, 58, 108]. Improving the prediction accuracy of thermal fluid behavior by CFD simulation is also essential. However, to evaluate the soundness of RPV in PTS, it is necessary to investigate until fracture analysis, and many such studies are being conducted (e.g., [10, 11, 123–132]).

Many analyses adopting the two-fluid model have been reported in the two-phase flow analysis of steam–water conditions. Therefore, the mesh resolution is 1 M cells (or elements) or less, the computational cost is kept low, and the turbulence model mostly uses RANS. Most studies used the interfacial momentum transport model AIAD for CFX and LIM for NEPTUNE_CFD. For LIM, the SSG model, suitable for highly anisotropic turbulence, was superior to the standard model [107], and many analyses using the SSG model have been reported. The DCC models of Hughes and Duffey [102], Coste [110], Banerjee et al. [122], Lakehal et al. [111, 121], and Magnaudet and Calmet [112] were considered as heat and mass transport models. Although it depends on the scope of the application [103], such as the turbulence intensity, the analyses validating the DCC model of Lakehal et al. [111] have been studied the most, and it has been shown that good accuracy can be obtained using it in combination with the SSG model [105]. The next most studied model is the DCC model based on the surface renewal theory proposed by Hughes and Duffey (HD model) [102], which is used in various codes. NEPTUNE_CFD showed excellent computational accuracy due to the accumulation of numerous validations. The analyses using the interface capturing/tracking method and the LES are challenging themes. Fundamental research focusing on single and combined phenomena should continue accumulating. Codes that can be applied not only to PTS phenomena but also to a wide range of reactor thermodynamic phenomena are required, and further research is expected.

Furthermore, we did not mention a small apparatus for benchmarking CFD codes for PTS here. However, in recent years, more detailed three-dimensional data obtained in the experimental apparatus for the CL mixing CFD-UQ benchmark [28] organized in the framework of the Working Group on Accident Management and Analysis (WGAMA) of the OECD/NEA and validating CFD simulation has been performed [133–140]. This experiment was performed under normal temperatures and pressures. However, the detailed data acquisition from high-temperature and high-pressure experiments for CFD simulation is still only halfway, and improvements in experimental measurement technology should expand the possibilities for CFD simulation validation.

Abbreviations

AIAD:	Algebraic interfacial area density
AFM:	Algebraic flux model
BMFT:	Bundesministerium für Forschung und Technologie
BSL:	Baseline
CATHARE:	Code for Analysis of Thermalhydraulics during an Accident of Reactor and Safety Evaluation
CFD:	Computational fluid dynamics
CI:	Core inlet
CL:	Cold leg
COSI:	Condensation due to safety injection
CRP:	Cooperative Research Project
DC:	Downcomer
DCC:	Direct contact condensation
ECC:	Emergency core cooling
ECCS:	Emergency core cooling system
EDF:	Électricité De France
GGDH:	General gradient diffusion hypothesis
HD:	Hughes and Duffey
HZDR:	Helmholtz–Zentrum Dresden–Rossendorf
IAEA:	International Atomic Energy Agency
IAPWS:	International Association for the Properties of Water and Steam
JAEA:	Japan Atomic Energy Agency
KKG:	Kernkraftwerk Gösgen
KWU:	Kraftwerk Union
LBLOCA:	Large break loss-of-coolant accident
LDV:	Laser Doppler velocimeter
LES:	Large eddy simulation
LIM:	Large interface model
LOCA:	Loss-of-coolant accident
LS:	Level set
LSTF:	Large-scale test facility
MC:	Magnaudet and Calmet
MS:	Mixing scalar
MSLB:	Main steam line break
NEA:	Nuclear Energy Agency
NURESAFE:	Nuclear Reactor Safety Simulation Platform
NURESIM:	Nuclear Reactor Simulations
NURISP:	Nuclear Reactor Integrated Simulation Project
OECD:	Organisation for Economic Co-operation and Development
PKL2:	Primärkreislauf 2
PS:	Pump simulator
PTS:	Pressurized thermal shock
PWR:	Pressurized water reactor
QUICK:	Quadratic upstream interpolation for convective kinetics
RANS:	Reynolds averaged Navier-Stoke
RELAP:	Reactor Excursion and Leak Analysis Program
SSG:	Second-order Speziale–Sarkar–Gatski Reynolds stress model
RNG:	Renormalization group
ROCOM:	Rossendorf coolant mixing model
ROSA:	Rig of Safety Assessment
RPV:	Reactor pressure vessel
RSM:	Reynolds stress model
SBLOCA:	Small break loss-of-coolant accident
SGDH:	Simple gradient diffusion hypothesis
SGS:	Subgrid-scale
SLR:	Standard low-Reynolds
SO:	Stuck-open valve
SRS:	Scale-resolving simulation
SSSW:	Steady-state steam–water test
SST:	Shear stress transport
TAMU:	Texas A&M University
TOPFLOW:	Transient two-phase flow
TSW:	Transient steam–water
UPTF:	Upper plenum test facility
UQ:	Uncertainty quantification
URANS:	Unsteady Reynolds-averaged Navier-Stoke
VLES:	Very large eddy simulation
VOF:	Volume of fluid
WALE:	Wall-adapting local eddy viscosity
WGAMA:	Working Group on Accident Management and Analysis
WMS:	Wire-mesh sensors.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was performed under the contract research entrusted by the Secretariat of Nuclear Regulation Authority (NRA), Japan. The authors would like to thank Dr. T. Takeda and Mr. K. Hiyama of JAEA for providing the literature on LSTF and a window person on software for supercomputers, respectively, Dr. S. Abe of NRA for providing the knowledge on CFD simulation for PTS, CEA/DM2S/STMF/LMSF; Dr. M.G. Rodio of CEA for providing the study opportunity on TOPFLOW-PTS and NEPTUNE_CFD cord; Dr. S. You of ANSYS Japan K.K. for answering the questions on software; and Dr. T. Watanabe of former Fukui University for providing the literature and experiment information. One of the authors (T. Hibiki) appreciates the support through a Global STEM Professorship in Hong Kong and the Hong Kong Jockey Club.