Abstract
The external reactor vessel cooling provides a means to stabilize the molten core inside the reactor vessel during the severe accident in nuclear power plants. For a credible evaluation of the chances of success of in-vessel melt retention, one requires high-fidelity predictions on the natural circulation flow around the reactor vessel and, in particular, the corresponding critical heat flux (CHF) on a hemispherical lower head. In this study, the coolability limit of external reactor vessel cooling in the Korean APR1400 reactor was investigated numerically by improving a thermal-hydraulic system analysis code, MARS-KS1.5. Since such system analysis codes did not incorporate a CHF model for a downward-facing hemisphere, three CHF correlations were newly implemented in the MARS-KS code, two of them taking into account the effect of the mass flux. A transient analysis on the CHF and surface heat flux on the lower head was carried out while the surrounding cooling water was heated to create the two-phase natural circulation flow. The heat transfer mode on the heated surface was used as a measure of evaluating the success of the strategy, and the simulation results were compared with those obtained via the default CHF model in the code. The sensitivity analysis on the heat load from the core melt was also performed.
1. Introduction
In a nuclear power plant, a prolonged undercooling of the nuclear reactor may cause the rupture of fuel rods and substantial meltdown of the core. The degradation of the reactor core allows the fission products that had been confined in the fuel rods to be released into the primary circuit. If the integrity of the primary piping and the pressure-bearing containment building should be further compromised, the severe accident may pose undue risk to public health on account of the unmitigated release of radioactive materials into the environment. In the aftermath of Fukushima Daiichi accident in 2011, many nuclear regulatory authorities have strengthened requirements for severe accident management of modern nuclear power plants so that the radiological consequence of the hypothetical severe accident can be managed.
The in-vessel retention is a scheme for severe accident management to arrest the relocation of molten corium at the lower head of the reactor vessel by flooding the reactor cavity, submerging the reactor vessel [1, 2]. This strategy is intended to maintain the lower head integrity via cooling its external surface to cope with the thermal load created by the high-temperature melt on the inside. For reactors having thermal insulation around the reactor vessel, an upward two-phase natural circulation flow will be established through the annular channel between the reactor vessel and insulation as described in Figure 1. If successful, the in-vessel retention strategy prevents uncontrolled release of fission products into the containment as well as energetic ex-vessel reactions, such as steam explosion and corium-concrete interactions. The in-vessel retention strategy was first incorporated into the Loviisa VVER-440 reactors in Finland [3] and thereafter extensively investigated for implementation into US AP600 designed by Westinghouse [4, 5].
The coolability limit of external reactor vessel cooling is defined by the critical heat flux (CHF). That is, for successful melt retention, the heat flux imposed by the molten corium in the lower head should be less than the CHF on the exterior surface of the lower head, which is a downward-facing hemisphere, to avoid boiling crisis [4]. It should be noted that the CHF is greatly influenced by the surface geometry of the heated wall. On a downward-facing hemisphere, the vertical migration of vapor bubbles is interfered with by the heated wall, which makes harder for bubbles to escape. Hence, the CHF on a downward-facing hemisphere is lower than that on a vertical wall, especially, at local locations with a lower polar angle (0° indicates the bottom of the lower head).
This indicates that a credible evaluation of the cooling capacity of external reactor vessel cooling requires reliable CHF data or models for the downward-facing hemisphere. Up to date, a number of experimental studies have been conducted to measure the CHF on the lower head according to the polar angle. The experimental facilities are basically divided into three categories: full-scale facilities with the slice geometry of the reactor vessel [1, 5–7], subscale facilities with a three-dimensional simulation of the reactor vessel [8], and separate-effect facilities with a flat rectangular channel whose inclination angle can be altered [9–11].
Once the CHF on the lower head is estimated, a sufficient margin between the CHF and the wall-to-fluid heat flux on the exterior wall should be preserved for the success of in-vessel retention. Because of huge thermal loads from the molten core, nucleate boiling will occur on the lower head wall as soon as external reactor vessel cooling begins. It should be noted that the nucleate boiling heat transfer coefficient is supposed to vary in time depending on flow characteristics formed by ex-vessel natural circulation, whereas the transient behavior of the corium pool will be much slower. Thus, an accurate prediction of nucleate boiling on a downward-facing heated wall is considered another challenge for assessment of the performance achievable by external reactor vessel cooling, and several recent publications also addressed this issue [12–14].
In most cases, the numerical evaluation of in-vessel retention through external reactor vessel cooling was carried out by using integral codes for severe accident analyses. Tarabelli et al. applied the ASTEC code to simulate the in-vessel retention with the improved modeling for thermal behavior of the core melt in the lower head [15]. The in-vessel module to describe the corium behavior (DIVA) and the ex-vessel module to model the external reactor vessel cooling (CESAR) were coupled in this study. The ASTEC V.12 code was applied to simulation of in-vessel retention of VVER-440/V213 reactor, but the convective boundary conditions on the ex-vessel wall were implemented in a simple way of imposing constant heat transfer coefficient and coolant temperature.
Jin et al. investigated the in-vessel retention strategy of a three-loop, 5000 MWt pressurizer water reactor using the MELCOR code (The code version was not addressed.) [16]. In this study, attempts were made to perform a coupled analysis that considers both internal and external conditions of the reactor vessel by modeling the flow path of the naturally driven cooling water outside of the reactor vessel. Large fluctuations in the coolant flow rate were reported after the transition to the saturated state, but they were merely explained as a matter of the numerical treatment of the code, and their effect on the CHF could not be reflected because the selected CHF correlation (ULPU-II correlation [6]) was only expressed in terms of the polar angle on the lower head.
Lim et al. simulated the severe accident progress of a Korean high-power reactor, APR1400, initiated by a large-break loss-of-coolant accident (LBLOCA) using MELCOR 2.1 [17]. Regarding heat transfer at the lower head, a simplified pool boiling model was used for nucleate boiling of cooling water on the ex-vessel wall, and the CHF was calculated by the El-Genk and Guo correlation [18], which is expressed in terms of fluid properties and the polar angle at the lower head. In this study, the cooling water in the reactor cavity was simply assumed to be at the saturation temperature, and the flow rate of cooling water and its impact on the CHF were not predicted as the cooling channel was modeled using a single control volume.
The literature reviews indicated that previous numerical research on in-vessel retention through external reactor vessel cooling has placed more attention on modeling the behavior of the molten core debris, whereas it relied on nucleate boiling models developed on pool boiling conditions or CHF correlations expressed as a function of the polar angle only in simulating the ex-vessel cooling. However, it should be noted that both convective heat transfer coefficient and CHF are impacted significantly by flow conditions over the heated surface. A recent experimental study to measure the CHF for in-vessel retention also reported that the CHF could shift by up to 35% according to the mass flux under the same thermal conditions (see Figure 12 in [11]). Hence, for a realistic prediction of the margin to the CHF during external reactor vessel cooling, (1) a high-fidelity prediction of the natural circulation flow established around the reactor vessel is required, and (2) a CHF model and nucleate boiling model for the downward-facing surface that can consider the flow effect should be incorporated in the analysis.
In this study, we simulated external reactor vessel cooling using the thermal-hydraulic system analysis code, which is a reliable tool for analysis of two-phase natural circulation flows. The target reactor was APR1400, and we selected the MARS-KS code as an analysis code, which is used for regulatory assessment by Korean nuclear authority [19]. An input model of the MARS-KS code to predict the time-dependent natural circulation flow rate by external reactor vessel cooling had been developed in our prior study via comparison with calculation results of a computational fluid dynamics (CFD) code. In this study, we improved the MARS-KS code by implementing three selected CHF correlations for the downward-facing hemisphere so that we could achieve a more systematic evaluation of the limit of coolability of external reactor vessel cooling in APR1400. For a specified heat load by the molten core, we performed the transient analysis to investigate the variable CHF and the heat transfer mode on the exterior surface of the lower head.
2. External Reactor Vessel Cooling of APR1400
2.1. Ex-Vessel Cooling System
The in-vessel retention through external reactor vessel cooling was employed in APR1400, a Korean advanced light water reactor, as a part of severe accident management strategy, even though it was not considered as a licensing design basis [20]. The water source for external reactor vessel cooling is the in-containment refueling storage tank (IRWST). Because of the configuration of the reactor cavity and the IRWST, the cooling water is initially supplied by the cavity flooding system as illustrated in Figure 2. Once the core exit temperature exceeds 650°C (1200°F), operators are allowed to actuate a shutdown cooling pump for flooding the reactor cavity up to a designated level. After the operation of the shutdown cooling pump is terminated, a boric acid makeup pump (BAMP) can be manually brought into operation to compensate for the loss of the water inventory by boiling on the ex-vessel wall.
As illustrated in Figure 1, the reactor vessel is surrounded by the insulation to minimize the heat loss from the reactor vessel under normal operation and limit temperatures of the reactor cavity concrete as well as ex-core sensors during accident conditions. Once the reactor cavity is flooded, the natural circulation flow channel is established around the insulation. That is, the cooling water heated by the thermal load from the relocated molten core moves upward through the inner annulus channel between the reactor vessel and the insulation. The steam/water mixture formed at the lower head is allowed to escape the outer annular gap between the insulation and the reactor cavity wall through egress doors at the side wall of the insulation. Then, steam is vented out of the reactor cavity to the containment atmosphere, while the separated water recirculates through the outer annulus channel.
The cooling water is supposed to enter the inner annulus channel through ingress doors installed at the bottom of the insulation. Normally closed, ingress doors are mostly arranged around the outer circumference of the bottom of the insulation in APR1400 [21]. Both of ingress panels and egress panels are hinged to the insulation. Note that the insulation was designed in a way that the flow paths along the inner annulus channel are narrowed around the lower head, thereby increasing the flow velocity; this will help somewhat enhance the CHF on the lower head.
2.2. Problem Definition
During in-vessel retention of the core melt, the reactor vessel may fail by a creep rupture of the reactor vessel by persisted thermal and pressure loads [22], the ejection of penetration tubes of in-core instrumentation at the bottom head, or the focusing effect of the heat flux on the inner surface of the reactor vessel wall induced by the metal layer [23]. Among them, this study is aimed at the failure by the focusing effect; if a corium pool forms in the lower head, the pool may be stratified with a metal layer on the top of an oxide pool, and the side wall in contact with the metal layer is then subject to a highly elevated heat flux.
Compared to previous research, the numerical study herein eliminated many assumptions associated with flow conditions inside the reactor cavity, such as constant water temperature, or constant heat transfer coefficient, or indefinite water inventories. Simulation conditions of this study were basically the same as those presented in [21]. The initial state of the simulation was defined as the situation in which a steady-state configuration of the core melt pool was formed inside the lower head and the cooling water filled up the reactor cavity as described in Figure 1.
The initial water temperature was set to 48.9°C. The atmospheric pressure of the containment building, and thus the pressure applied on the water surface, was set to 250 kPa considering the elevation of pressure caused by a hypothetical release of the reactor coolant. In this simulation, it was assumed that the manual operation of the BAMP began as soon as the collapsed water level in the outer annulus channel was lower than the elevation 0.91 m below hot legs. Then, the BAMP delivered make-up water from IRWST at a designated flow rate of 38.6 m3/hr (170 GPM).
The controlling boundary condition for external reactor vessel cooling is the thermal load posed by the molten corium from inside the lower head. In this study, the heat flux profile calculated by the MAAP4 code for the core degradation scenario of APR1400 initiated by the large-break LOCA was implemented as the boundary condition for the inner wall of the lower head. The applied heat flux distribution according to the polar angle is shown in Figure 3. The heat flux increased with the polar angle, and the maximum value reached 1.27 MW/m2 over the region of the light metal layer, which rode on the top of the heat-generating oxide pool. The corresponding total thermal load was 23.3 MW.
3. Selection of CHF Correlations
In order to select the most appropriate CHF correlations for a new simulation approach using MARS-KS, we extensively reviewed past experimental studies in which the CHF for in-vessel retention had been measured. Our primary emphasis was on selecting and implementing CHF correlations that can predict the flow effect on the CHF. A considerable effect of the mass flux on the CHF has been already demonstrated in several experimental studies even for the downward-facing hemisphere geometry; it is of greater importance for high-power reactors, such as APR1400, since the margin to the CHF is expected to be small. From the literature review, we selected three CHF correlations to be implemented into the MARS-KS code and tested: ULPU-II correlation [4], Toshiba correlation [11], and SULTAN correlation [10]. The test sections of experimental facilities for selected test programs are shown in Figure 4.
(a) ULPU-II test
(b) Toshiba test
(c) SULTAN test
The ULPU facility was originally built to simulate AP600 geometries and ex-vessel cooling conditions and extended to investigate the performance of a higher-power reactor, AP1000. The ULPU facility employed a full-scale slice model of the cooling channel outside of the reactor vessel and hemispherical-slice heater blocks. A total of five test programs were performed varying the configuration of the flow circuit along the hemispherical head, which was controlled by the presence and the design of a baffle simulating the lower part of the insulation. Among them, the ULPU Configuration II introduced no baffle while ULPU Configurations III–V installed a flat or a curved baffle to induce a higher flow velocity by natural circulation. The lower envelope for test data from Configuration II, in kW/m2, was expressed as a function of the polar angle as (Let us call it the ULPU-II correlation in this paper):
Even though Equation (1) does not account for the flow effect on the CHF on a lower head, it was chosen because it is the most cited correlation for the evaluation of in-vessel retention and is likely to provide a somewhat conservative estimate.
Recently, researchers at Toshiba Energy Systems & Solutions Corporation carried out an experimental study to newly obtain CHF data for in-vessel retention equivalent to actual plant conditions in pressurized water reactors [11]. The Toshiba correlation was selected since its experimental work systematically elucidated the effect of the mass flux. To simulate the ex-vessel cooling channel, a rectangular test section with a cross-section of and a heated length of 600 mm was used. The heater block made of copper was installed on the upper side of the rectangular flow channel. The inclination angle of the test section was set to 50°, 70°, and 90° considering the lower and the upper bounds of positions of the light metal layer in the stratified molten corium. The flow conditions in the test section were controlled by forced convection, and a total of 36 test cases were obtained to figure out the parametric effect of pressure, mass flux, and thermodynamic quality.
The CHF in kW/m2 was derived as follows (Let us call it the Toshiba correlation in this paper): which is valid for pressure ranging from 0.1 to 0.6 MPa, mass flux from 430 to 1900 kg/m2s, thermodynamic quality from to , and polar angle from 50 to 90°.
The SULTAN program run by CEA was also intended to increase CHF data in a wide range of thermal-hydraulic and geometric parameters [9, 10]. The heated section is a flat plate, electrically heated, of 4 m in length and 15 cm in width in a rectangular channel. The channel width could be varied from 3 to 15 cm, and the test section could be inclined from the horizontal to the vertical position. Two-phase flow parameters were also measured in the test, including the pressure drop, void fraction distributions, and temperature distributions. A total of 191 CHF data were obtained from the SULTAN facility. The SULTAN correlation was proposed in terms of the pressure, mass flux, thermodynamic quality, gap width, and inclination angle to estimate the CHF in MW/m2 as follows:
Details of coefficients – in Equation (3) are referred to [10]. The valid range of the SULTAN correlation lies within 10–90° for the inclination angle, 0.1–0.5 MPa for the pressure, 10–5000 kg/m2s for the mass flux, and 0–50 K for the inlet subcooling. The SULTAN correlation was chosen since it had merits of not only considering various flow and geometric parameters but also covering the wide range of flow conditions, which helps robust calculation of the CHF when implemented into an analysis code.
In open literatures about the above tests, either a full set of test conditions (in terms of , , , ) was not provided, or only a limited portion of measured CHF data was available. Hence, merely an indirect comparison among experimental results was possible by plotting the calculated CHF from proposed correlations. Figure 5 presents the predicted CHF versus the coolant mass flux on the vertical position of the lower head () at 0.3 MPa. The thermodynamic quality was assumed to 0, and the gap width was set to 10 cm. Since the ULPU-II correlation considers only the variation of the CHF with the polar angle, the calculated CHF was irrelevant to the coolant mass flux. Thus, we cannot achieve a realistic simulation of thermal margins by relying on the ULPU-II correlation in an aspect that the effect of the natural circulation flow established around the lower head is ignored. Both the Toshiba correlation and SULTAN correlation predicted the increase of the CHF with the coolant mass flux. In the saturated liquid state, the CHF predicted by the SULTAN correlation was lower than that predicted by the Toshiba correlation by about 30%.
4. MARS-KS Model
4.1. Field Equations and Closure Models
The MARS-KS code is a nuclear safety analysis code used for the thermal-hydraulic analysis of all transients and postulated accidents in light-water reactor systems [19]. In MARS-KS, the two-phase flow is described by a one-dimensional transient two-fluid model. The continuity equations, momentum equations, and energy equations are solved for liquid and vapor phases, respectively, as follows:
For closure of field equations, constitutive models and correlations are incorporated for processes such as interphase friction, interphase heat transfer, wall friction, and wall heat transfer. Among them, the CHF model is used for the identification of not only the two-phase flow regime but the heat transfer mode on a solid surface on which a convective boundary condition is imposed. The default CHF model of the MARS-KS code, which is intended for an internal flow in circular pipes, adopted the 1986 AECL-UO Critical Heat Flux Lookup Table method [24]. The lookup table was formulated from test data normalized to a tube of 0.008 m in inner diameter, and eight multipliers are introduced to consider the effects of the tube diameter, rod bundles, flow directions, axial power shape, and so on.
For the default geometry, the heat transfer coefficient of the single-phase liquid is determined as the maximum of two values calculated for the turbulent forced convection (Dittus-Boelter equation) and the free convection (Churchill-Chu correlation for a vertical cell, McAdams correlation for a horizontal cell), respectively [25].
While the CHF calculated by the implemented correlations defines the coolability limit on the heated wall, the nucleate boiling heat transfer coefficient determines the actual wall-to-fluid heat flux on the outer surface of the wall in the “pre-CHF” regime. Influenced by flow conditions around the lower head, the nucleate boiling heat transfer coefficient is also used to solve the unsteady conduction equation for heat structures to determine the surface temperature in time as well as the temperature profile of the wall. For most hydraulic geometries including the downward-facing hemisphere, the nucleate boiling heat transfer is predicted by the Chen correlation [26].
4.2. Implementation of New CHF Correlations
Three selected CHF correlations described in Section 3 were implemented into MARS-KS version 1.5 by modifying its source code. The subroutines to calculate the CHF based on each correlation were added so that they could be called when the heat transfer advancement on the downward-facing hemisphere is required. In the MARS-KS code, users are supposed to specify the geometry type in the input deck when they pose a convective boundary condition for the heat structure. Then, the corresponding heat transfer package is called to classify the heat transfer mode and calculate the wall-to-fluid heat transfer rate.
The convection boundary types that can be selected in the MARS-KS code are listed in part in Table 1. In this study, convection boundary types were extended to include the description of the downward-facing curved surface through indices 190–192. That is, the geometry type number 190 was designated to call the ULPU-II correlation, 191 the Toshiba correlation, and 192 the SULTAN correlation in the modified MARS-KS code. The source code was also modified in a way that the polar angle of a heat structure and the gap width, just in case the SULTAN correlation is called, are read from the input cards once the geometry type index belongs to 190–192. We confirmed that these modifications caused no interference with embedded modules for heat transfer calculations.
Note that valid flow ranges for the Toshiba correlation are somewhat limited since the test matrix had been based on the natural circulation flow rate estimated under the balanced conditions. In addition, the correlation focused only on high polar angle ranges between 50 and 90° where the focusing effect by the light metal layer is concerned. Our preliminary calculations revealed that the Toshiba correlation must not be extrapolated out of the valid ranges; the estimated CHF was very low under low mass flux conditions at low thermodynamic quality of cooling water. Hence, the Toshiba correlation was coded in a way that, if a flow parameter deviated from the valid range, the bounding value of the valid range of the parameter replaced the actual value in computing the CHF to insure the stability of the calculation.
4.3. Input Model
The MARS-KS nodalization of the natural circulation flow channel during external reactor vessel cooling of APR1400 is illustrated in Figure 6. The flow circuit was modelled by a total of 56 nodes and 57 junctions. The ingress openings at the bottom of the insulation were represented by the single-junction component C205. The annular gap between the lower head and the insulation was modelled by three annulus components C220, C230, and C240, comprising a total of 10 nodes. These volumes represent the effective heating section by the core melt hypothetically filling the reactor cavity. The heat load to the coolant was simulated by connecting the heat structures to C220–C240, i.e., S220–S240. The total heat rate into each node was obtained via the integration of the heat flux profile in Figure 3 over the corresponding range of the polar angle.
C345 was used to model the egress openings through which the steam/water mixture is vented to the outer annular channel between the insulation and the reactor cavity wall. The downward recirculation channel of water was modelled by annulus components C410 and C420, while steam is supposed to escape through C510. The pressure boundary condition was implemented by the time-dependent volumes C330 and C520, which simulate the atmosphere of the containment. The coolant makeup by the BAMP was modelled by the time-dependent junction C115 connected to the third hydrodynamic volume for the lower reactor cavity. More details of the nodalization for the MARS-KS analysis are referred to [21].
Due to the nature of one-dimensional field equations, such system-scale analysis code as MARS-KS has some limitations in fully describing local flow characteristics induced by complicated geometries. Especially, in case of simulating free convection problems with this kind of codes, the calculated natural circulation flow rate is highly dependent on user-supplied loss coefficients in the input model because the driving force for free convective flows is normally very small. However, we do not have a standard method to determine hydraulic resistances of flow through such large-scale annular channel with varying cross-sectional area and localized openings as illustrated in Figure 1. To minimize the uncertainty of flow rate calculations originated from the ambiguity in local loss coefficients, we carried out in our recent study a three-dimensional CFD (computational fluid dynamics) analysis to estimate the natural circulation flow rate by external reactor vessel cooling. Details of the CFD analysis are explained in [21].
From the CFD simulation, we estimated the flow rate and pressure distribution of a liquid-state coolant just before significant bubble generation as shown in Figure 7. They were used to obtain the loss coefficient at three flow junctions where substantial minor losses are expected: the ingress openings (C205), the egress openings (C345), and the downward channel between the lower part of the insulation and reactor cavity wall (C420-3). The determined loss coefficients rendered the flow rate calculated by the MARS-KS code identical to that predicted by the CFD analysis when the inlet temperature of cooling water into the inner annulus channel reached a designated value. The loss coefficients obtained by the above CFD-aided estimation were given in the input deck of the MARS-KS code.
(a) Velocity magnitude
(b) Absolute pressure
5. Results and Discussions
5.1. Coolability Limit and Heat Transfer Mode
Using the modified MARS-KS code, transient analyses were conducted on the case in which the thermal load from the molten corium totaling 23.3 MW was applied to the cooling water filling up to the bottom of the cold leg. The end time of the problem was set to 10,000 seconds. As investigated in the previous analysis, the cooling water reached the saturation temperature, and a quasi-steady two-phase natural circulation flow was developed long before 10,000 seconds elapsed.
Figures 8 and 9 presents the calculation result of the coolant temperature around the lower head and the natural circulation flow rate through the ingress openings (C205), respectively, from the input model reflecting the obtained loss coefficients. As high heat loads were imposed to the cooling channel, the temperature of the cooling water continued to rise and approached the saturation value. While the cooling water was at the subcooled state, the flow rate was not more than 300 kg/s. As the subcooling degree was reduced below a level enough to cause significant vapor generation on the lower head, the natural circulation flow rate increased drastically. Note that, during the transition into the two-phase natural circulation flow regime proceeded, significant oscillation of the flow rate occurred by the two-phase flow instability. Then, the natural circulation flow was stabilized as the cooling water around the lower head was sustained at the saturated state. Conventionally, the chances of success of external reactor vessel cooling had been investigated at this two-phase regime (phase C in Figure 9).
Figure 10 shows the CHF and surface heat flux on the exterior surface of the lower head predicted by the original MARS-KS code with the default CHF lookup table. The simulation results in Figure 10 are those of the uppermost heat structure of the effective heating section, S240–2, where the heat flux from the molten corium was the highest and the temperature (or quality) of the cooling water was maximum. It was shown in Figure 10 that the surface heat flux was lower than the CHF by more than 1 MW/m2 at the subcooled state. However, with an oscillation of the flow rate in the transition regime, both the surface heat flux and the CHF also fluctuated substantially. In the meantime, the surface heat flux temporarily exceeded the CHF, causing the transition of the heat transfer mode into the transition boiling regime. By the time the flow oscillation by free convection was terminated, the gap between the surface heat flux and the CHF widened again. It seemed as if the strategy of external reactor vessel cooling could effectively prevent the thermal failure of the reactor vessel.
However, the above analysis was attributable to the fact that the CHF calculated by the default model of the MARS-KS code was the value in the vertical tube. That is, the default model overestimated the CHF on the downward-facing hemisphere. Hence, the above analysis relying on the default model was not realistic, nor conservative in a view of nuclear safety. This demonstrated the reason why the empirical correlations based on the experimental data for the downward-facing hemisphere must be applied for the thermal evaluation of external reactor vessel cooling.
The modified MARS-KS1.5 code was run with the identical input model described in Section 4.3 to evaluate the limit of coolability of external reactor vessel cooling according to the selected CHF correlation. Figures 11 and 12 show the calculation results of the CHF and surface heat flux on the heating surface of S240–2, influenced by the evolution of the natural circulation flow rate, when the Toshiba correlation and the SULTAN correlations were applied, respectively. The time-independent CHF predicted by the ULPU-II correlation was also drawn in Figure 11.
Note that, in the early phase of external reactor vessel cooling, the predicted CHF by the SULTAN seems larger than that by the Toshiba correlation. This is associated with the way each correlation was implemented in the MARS-KS code. The subcooling degree of the cooling water in the reactor cavity was initially as high as 85 K, and out of the valid range of two correlations. While the SULTAN correlation was extrapolated for the subcooling degree over 50 K, the lower limit of the thermodynamic quality () was imposed for the Toshiba correlation as described in Section 4.2. This was the reason why the calculated CHF in Figure 11 was relatively flat until 1,500 s.
The simulation results revealed that, in any case applying a CHF correlation for the downward-facing hemisphere, the surface heat flux exceeded the CHF repeatedly due to flow oscillations in the early stage of the transition regime alike the simulation result with the default CHF model. Figure 13 presents the transition of the heat transfer mode on S240–2 per the selected CHF correlation. While the cooling water was in the subcooled state, the heat removal from the lower head occurred in the nucleate boiling regime. Then, the repeated transitions between nucleate boiling and transition boiling were observed in the flow oscillation phase.
The rate of convective heat transfer is retarded in the transition boiling regime. Hence, the longer that the lower head wall underwent transition boiling, the wall temperature increased rapidly to lead to the transition to film boiling, eventually, in all cases. The flow chart to determine the post-CHF heat transfer regimes is presented in Figure 14 [25]. The MARS-KS code regards directly that post-CHF heat transfer occurs on the heated surface if the surface temperature exceeds the saturation temperature plus 600 K and no longer requires the CHF for comparison with the surface heat flux. This is the reason why the CHF was not calculated after about 2,130 s in simulations of Figures 11 and 12.
Note that the complete transition to film boiling occurred in less than a minute after the first experience of transition boiling. As the heated surface experiences the departure from nucleate boiling (DNB), the MARS-KS code determines the post-CHF heat transfer mode based on the relative magnitude of the heat flux estimated from the transition boiling correlation and the film boiling correlation as shown in Figure 14. For film boiling, the MARS-KS code adopted the correlation developed by Bromley [27], in which the heat flux increases with the difference between the wall temperature and the saturation temperature. Due to the huge heat source to S240–2, the rate of the wall temperature rise was so high that the heat flux calculated by the film boiling model increased very fast, and transition boiling sustained temporarily.
This is contrary to the simulation result by the default CHF model in which post-CHF heat transfer did not last long, and the heat transfer mode returned to nucleate boiling soon. When the SULTAN correlation was applied, the surface temperature of the lower head sustained below 446 K in nucleate boiling turned out to rise over 1,500 K in 250 seconds after transition boiling appeared in the flow oscillation phase as shown in Figure 13(c). Since the MARS-KS code is incapable of simulating the melting of heat structure materials, the transition of the heat transfer mode to film boiling can be interpreted that the thermal load from the molten core exceeds the limit of coolability of external reactor vessel cooling, thereby resulting in thermal failure of the lower head. The above results demonstrated that not only the CHF itself but also the fluctuation of nucleate boiling heat transfer in the flow oscillation phase, as shown in Figures 11 and 12, have a critical impact on the thermal margin to the lower head failure.
To understand the frequent transition of boiling regimes in the flow oscillation phase, it is required to review the logic of selecting the heat transfer mode and the corresponding wall-to-fluid heat transfer models in the MARS-KS code. In the MARS-KS code, the heat structure is considered to experience the DNB if the wall temperature is greater than the coolant saturation temperature by over 100 K or the heat flux is larger than the CHF as shown in Figure 14. If either of these two requirements is not met, then the heat transfer mode is determined to nucleate boiling. The selection of the heat transfer mode is performed at each time step.
In the subcooled natural circulation regime (phase A in Figure 9), the temperature difference between the wall and the fluid at C240–2 was maintained at more or less 50 K. In the early stage of the transition regime (phase B), the surface heat flux continued to change so rapidly due to natural circulation flow oscillations that the temporary excess over the CHF was repeated regardless of the applied CHF model. In the case of applying the default CHF model, the transition boiling regime lasted only short since the thermal margin to the CHF was large except for the early stage of phase B. Once again, this resulted from using the CHF database inappropriate for the lower head, overestimating the limit of coolability.
On the contrary, it was revealed from applying CHF correlations for the downward-facing hemisphere that the margin between the surface heat flux and the CHF was actually much smaller, and thus, the duration of post-CHF heat transfer was prolonged in the flow oscillation phase. Hence, the surface heat flux dropped as the heat structure began to undergo post-CHF boiling, in which the heat transfer coefficient is much less than that in nucleate boiling. The MARS-KS code predicted as if the heat flux to the coolant increased gradually again, but this arose simply from a continuous increase in the wall temperature without considering changes of state of the wall material; in a real situation, such an overheating would cause the lower head of the reactor vessel to fail.
5.2. Oscillation of the Surface Heat Flux
As a result, the above unexpected finding about the ex-vessel cooling performance originated from variability in the surface heat flux led by the natural circulation flow instability. Since the heat transfer model of nucleate boiling was not altered for the geometry types 190–192 in the modified MARS-KS, the heat transfer coefficient in the nucleate boiling regime was calculated by the Chen correlation [25, 26]. Figure 15 presents the variation of the cell-centered flow rate (mean value of inlet and outlet junction mass flows) at C240–2 and the calculated convective heat transfer coefficient. When the time interval for displaying the calculation results was limited to the initial 5 seconds in the flow oscillation phase, it is possible to observe in Figure 15(d) substantial oscillations of the convective heat transfer coefficient with time.
(a)
(b)
(c)
(d)
To figure out the oscillation of the surface heat flux in the early stage of phase B, the parametric effect of the flow rate on the heat transfer coefficient in forced-convective nucleate boiling was analyzed. The Chen correlation for saturated and subcooled nucleate boiling consists of a macroscopic convection term and a microscopic boiling term as where the forced convection component was calculated from the Dittus-Boelter equation, and the nucleate boiling component is based on the Forster-Zuber correlation as follows: where and . The factor in Equation (7) accounts for the enhanced flow and turbulence due to the presence of vapor. It is expressed as where
The suppression factor is a function of the total Reynolds number as where
We chose two temporal points Figure 15(c) to look to the variation of the heat transfer coefficient with flow conditions: one is the instant on the brink of flow oscillations, and the other is the instant at which the DNB occurred for the first time. These two temporal points were deviated by 1.07 seconds. The wall superheat changed little between two points (~52 K) at 3 bar, and the flow was in the low quality region. The cell-centered mass flux dropped rapidly from 89.9 to 25.6 kg/m2s, and thus, the flow quality at the latter point was about four times larger than at the former point. Under these conditions, the total heat transfer coefficient was controlled by the microscopic boiling term. It was revealed from the analysis that the suppression factor soared at the latter point by 85%, and the total heat transfer coefficient increased by 73%. This indicates that the enhanced heat transfer by two-phase nucleate boiling was crucial. In nucleate boiling at low quality, a decrease in the flow rate induced vigorous vapor generation to render the contribution of the microscopic boiling mechanism predominant.
Note that this study was motivated by a lack of considering the change in the CHF on the lower head by the natural circulation flow rate. However, the above numerical analysis revealed that the oscillation of the wall-to-fluid heat flux arising from the natural circulation flow instability could have a greater impact on coolability limits of external reactor vessel cooling. It was found from this study that an abrupt drop in the flow rate in the low quality, flow oscillating regime could lead to a drastic increase of the surface heat flux on the lower head, exceeding temporarily the limit of coolability of ex-vessel cooling strategy. We arrived at a conclusion that it is never sufficient to evaluate the thermal margin of external reactor vessel cooling based on flow conditions after a stable two-phase natural circulation flow is established as in phase C in Figure 9. The above results imply that it is also essential to evaluate whether the oscillations of the flow rate and the surface heat flux during the transition to the two-phase natural convection regime (phase B) lead to DNB or not.
In this study, we applied three CHF correlations for the downward-facing hemisphere as summarized in Section 3. According to our extensive preliminary calculations, the SULTAN correlation turned out to provide a robust calculation because it is valid over a wide range of flow conditions. If we also take into account slightly conservative estimates of the CHF as shown in Figure 5, the use of the SULTAN correlation (geometry type 192) is recommended for the simulation of external reactor vessel cooling with the modified MARS-KS code if a user has to pick just one.
5.3. Sensitivity Analysis
In the aforementioned problem, the total thermal load from the core melt was set to 23.3 MW. This was a representative value of the decay heat used in numerical studies on external reactor vessel cooling of APR1400 [28]. We performed a sensitivity analysis on the magnitude of the thermal load to investigate its effect on the unsteady natural circulation flow rate and consequential heat transfer mode on the lower head. This helps to estimate the upper limit of the corium heat load under which the boiling crisis does not occur on the lower head in spite of the instantaneous surge in the surface heat flux in the flow oscillation phase, and thus, the external reactor vessel wall can be cooled down.
In the sensitivity test, the total thermal load was reduced by 10% and 20%, respectively. That is, we assumed two scenarios with a heat load of 21.0 MW (case A) and 18.6 MW (case B). When the thermal load was raised over the reference value, the increase in the wall-to-fluid heat flux resulting from a larger heat source was predominant over the effect of the increased natural circulation flow rate, causing the boiling crisis on the lower head. The local heat flux versus the polar angle on the lower head was assumed to reduce from the reference profile shown in Figure 3 by the same ratio. The applied heat flux profiles in the sensitivity test are presented in Figure 16.
Figure 17 shows the unsteady cell-centered mass flow rate at C240–2 calculated in cases A and B. Like the simulation result of the reference problem, a substantial oscillation of the natural circulation flow rate was predicted during the transition into two-phase flow conditions, and a stable convective flow rate was established in the two-phase regime. A consistent decrease in the two-phase natural circulation flow rate was expected with a reduction in the thermal load from the core melt.
The CHF calculated from the SULTAN correlation and the surface heat flux in cases A and B were plotted in Figure 18. Both cases underwent abrupt changes in the surface heat flux in the flow oscillation phase, but the surface heat flux did not exceed the CHF. Accordingly, it was predicted that the lower head wall, including the heat structure connected to C240–2 as shown in Figure 19, could be cooled in nucleate boiling for the entire problem time without reaching DNB conditions. This sensitivity test implied that the thermal load set in the reference problem (23.3 MW) was close to an upper limit to which the coolability of external reactor vessel cooling in APR1400 could be attained to achieve in-vessel melt retention.
(a)
(b)
6. Conclusions
For a systematic evaluation of the coolability limits on the lower head achievable via external reactor vessel cooling, a Korean thermal-hydraulic system analysis code was improved by implementing three CHF correlations for the downward-facing hemisphere. Based on the developed input model to better predict the natural circulation flow rate established around the reactor vessel, the improved code was applied to a transient simulation of external reactor vessel cooling in APR1400.
It was found that applying the default CHF model of the thermal-hydraulic system analysis code can mislead the margins to thermal failure of the reactor vessel during external reactor vessel cooling since it overestimates the CHF on the lower head. Our transient analysis also revealed that it is crucial to explore whether the instantaneous surge in the surface heat flux caused by the oscillating flow rate leads to DNB during the transition to the two-phase natural convection regime.
The improved MARS-KS model from this study will be used to develop a new integrated analysis framework for in-vessel melt retention. That is, the MARS-KS model for ex-vessel cooling will be coupled with the smoothed-particle hydrodynamics (SPH) code developed by Seoul National University, named SOPHIA [29], which can simulate the behavior of the molten core inside the lower head based on the first principle. In the new analysis tool, the heat flux data from the molten core to the coolant calculated by the SOPHIA code are transferred to the MARS-KS code; then, the MARS-KS code calculates the exterior surface temperatures of the lower head and sends them to the SOPHIA code for advancement of heat flux profiles. This enables to simulate bidirectional interactions between the internal and external conditions of the reactor vessel wall during in-vessel melt retention. Our recent progress about codes coupling via the socket communication method and preliminary full-scale analysis of the in-vessel retention in APR1400 are described in [30].
Nomenclature
| : | Cross-sectional area (m2) |
| : | Body force (m s-2) |
| : | Specific heat at constant pressure (J kg-1K-1) |
| : | (Hydraulic) diameter (m) |
| DISS: | Energy dissipation function (W m-3) |
| : | Gap (m) |
| : | factor in the Chen correlation |
| : | Virtual mass force (kg m-2s-2) |
| : | Interfacial drag force (kg m-2s-2) |
| : | Wall frictional force (kg m-2s-2) |
| : | Mass flux (kg m-2-1) |
| : | Latent heat (J kg-1) |
| : | Macroscopic convection heat transfer coefficient (W m-2K-1) |
| : | Microscopic boiling heat transfer coefficient (W m-2K-1) |
| : | Thermal conductivity (W m-1K-1) |
| : | Pressure (Pa) |
| Pr: | Prandtl number, |
| : | Heat flux (W m-2) |
| : | Phasic wall heat transfer rate per unit volume (W m-3) |
| : | Interface heat transfer rate per unit volume (W m-3) |
| Re: | Reynolds number, |
| Retp: | Total Reynolds number |
| : | Suppression factor in the Chen correlation |
| : | Time (s) |
| : | Temperature (K) |
| : | Specific internal energy (J kg-1) |
| : | Velocity (m s-1) |
| : | Thermodynamic quality. |
| : | Void fraction |
| : | Volumetric mass exchange rate (kg m-3s-1) |
| : | Polar angle on the lower head (degree) |
| : | Viscosity (kg m-1s-1) |
| : | Density (kg m-3). |
| CHF: | Critical heat flux |
| : | Liquid phase |
| FB: | Film boiling |
| : | Gas phase |
| : | Phase |
| sat: | Saturated state |
| TB: | Transition boiling |
| spt: | Saturation based on total pressure |
| w: | Wall. |
Data Availability
Data are available on request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the 2022 scientific promotion program funded by Jeju National University.