Thermal-hydraulic analysis of a VVER-1000 core in MSLB conditions

. The objective of this paper is to analyze the ability of a VVER- 1000 core and its control system to cope with a hypothetical main steam line break (MSLB) accident in case of multiple equipment failures. The study involves the use of advanced 3D core calculation models benchmarked and validated for reactivity accidents in preceding studies. A MSLB core boundary condition problem is solved on a coarse (nodal) mesh with the coupled COBAYA/CTF neutronic/thermal hydraulic codes. The core thermal-hydraulic boundary conditions are obtained from a preceding full-plant MSLB simulation. The assessment of the core safety parameters is supplemented by a fine-mesh (sub-channel) thermal-hydraulic analysis of the hottest assemblies with the CTF code using information from the 3D nodal COBAYA/CTF calculations. Thirteen variants of a pessimistic MSLB scenario are considered, each of them assuming a number of equipment failures aggravated by eight control rods stuck out of the core after scram at different locations in the overcooled sector. The results (within the limitations of the adopted modeling assumptions) show that the core safety parameters do not exceed the safety limits in the simulated aggravated reactivity accidents. Further analyses expected include a coupled pin-by-pin neutronic/thermal-hydraulic calculation, at least for a mini-core of seven assemblies around the hottest one, as well as a dynamic gap conductance coefficient, a DNBR from VVER-specific CHF skeleton tables in CTF, and uncertainty quantification.


Introduction
A major concern in a main steam line break (MSLB) reactivity accident is the risk of core overheating. In the computational analysis of such accidents the safety parameters of particular interest are the fuel centerline temperature, the departure from nucleate boiling ratio (DNBR) and the fuel rod cladding temperature. The fuel temperature should be kept well below the UO 2 melting temperature which can significantly vary depending on the fuel burnup. As a VVER-1000 core can contain once, twice, three or four times burnt fuel assemblies, it is of practical interest to analyze the consequences of such an asymmetric reactivity accident at the end of core life and for various combinations of equipment failure. This paper presents results from the analysis of thirteen variants of a hypothetical MSLB scenario, each of them assuming a number of equipment failures plus eight control rods stuck out of the core after scram at different radial locations in the overcooled sector.
The main objective is to demonstrate the use of state-of-the-art 3D core calculation models which have been benchmarked and partly validated for VVER-1000 reactivity accidents in preceding studies. A specific objective is to make a step in the analysis of the ability of a VVER-1000 core and its control system to cope with such an accident without exceeding the safety limits. Uncertainty analysis is beyond the scope of this work.

Accident scenario
The considered accident scenario is based on aggravated variants of the pessimistic scenario in the OECD/NEA VVER-1000 MSLB benchmark (V1000CT-2) [1]. The task is to solve a MSLB core boundary condition problem using coupled 3D neutronic/thermal hydraulic codes, given the core thermal-hydraulic (TH) boundary conditions as obtained from a full plant simulation. The core boundary conditions (BCs) are taken to be as specified in the V1000CT-2 benchmark [1,2]. The plant transient is initiated at hot full power by a large guillotine type break of steam line #4 outside the containment, upstream of the steam intercept valve. The reference core is a real Kozloduy-6, Cycle 8 three-year batch core at 270.4 EFPD [1]. The core contains once, twice and three times burnt UO 2 fuel of initial enrichment 4. The steam generator in the faulted loop continues uncontrolled cooling till the complete evaporation of the secondary water. A cooler sector is formed at the core inlet, with overcooling of up to 80°C. During the transient eight peripheral control rod clusters (CR) are assumed to remain stuck out of the core after scram, all of them in the overcooled sector. Thirteen cases with different radial combinations of the stuck rods are to be analyzed to assess the values of the core safety parameters.

Full-core nodal calculation
Full-core coarse-mesh (nodal) simulations for each CR configuration were carried out with the coupled COBAYA4/CTF neutronic/thermal-hydraulic codes.
COBAYA4 [3,4] is a 3D multi-scale core physics code using transport-corrected multigroup diffusion approximation. It is developed by the Universidad Politecnica de Madrid and benchmarked for VVER-1000 calculations [5,6] in the frame of the EU NURISP [7] and NURESAFE [8] projects. At the nodal level the analytical coarse-mesh finitedifference (CMFD) method [9] is used. The code has radial mesh refinement capability.
COBRA-TF (CTF) [10,11] is a recent version of the COBRA-TF thermal-hydraulic (TH) code which uses a two-fluid, three-field modeling approach and has sub-channel capabilities.
The COBAYA4/CTF coupling method [12] for VVER-1000 is based on the MED Coupling libraries in the Salome platform [13]. The coupling and the coupled models have been tested for VVER MSLB in preceding studies [5,14]. The modeling assumptions in the coupled COBAYA4/CTF VVER-1000 calculation models are briefly summarized below: Coarse-mesh COBAYA 3D neutron kinetics with: -30 axial nodes in the heated region; -2 nodes in each axial reflector; -6 triangles per hexagon.
Use of a realistic VVER-1000 multi-parameter two-group cross-section library [15] for reactivity accident analysis which has been generated with the APOLLO2 code [16] and validated in the frame of the NURESAFE project [8]; Use of time-dependent assembly-by-assembly MSLB thermal-hydraulic core boundary conditions (inlet temperatures, inlet mass flow rates, and outlet pressures) obtained from a full plant simulation involving a quasi-3D reactor pressure vessel TH model [17]; Coarse-mesh CTF thermal-hydraulic model with one channel per assembly and 30 axial nodes in the heated part; Fuel model with 9 radial rings in the fuel, one for the gas gap and one for the cladding; Temperature-dependent fuel and cladding thermal-physical properties [18]; The spacer grids are not explicitly modeled and are taken into account by the vertical pressure loss coefficients; The gas gap conductance coefficient is taken constant, equal to 3070 W/m 2 K, as estimated at average core burnup of 26.6 MWd/kgHM [1]; Chen's model of nucleate boiling [19] and the W-3 critical heat flux (CHF) correlation [20] with non-uniform power distribution.

Sub-channel assembly calculations
Two variants of the transient having the most risky values of the core safety parameters were selected for sub-channel TH calculations of the hottest assemblies. Such simulations were carried out for the hottest assemblies and assemblies next to them so that the analysis includes once, twice or three times burnt fuel. The main modeling assumptions were as follows: Assembly thermal-hydraulic problems were solved with inlet/output TH BCs from the plant system simulation, and assembly powers and axial power profiles as obtained from the full-core nodal simulation with COBAYA4/CTF; Radial pin-power distribution taken such as at hot full power for the hottest assemblies, and with artificially imposed 5% radial tilt for the considered adjacent assembly of higher burnup; Coolant-centered radial spatial mesh with 660 sub-channels per assembly; Axial mesh with 30 nodes in the heated region; Fuel model with 9 radial rings in the fuel pellet, 1 for the gap and 1 for the cladding. The central hole is taken into account, and conduction in radial and axial direction is considered; The bypass of 2.2% through the control rod guide tubes in the un-rodded assemblies is not explicitly modeled and is taken into account by decreasing the active coolant flow; The spacer grids are not explicitly modeled and are taken into account by the vertical pressure loss coefficients; User defined coolant mixing coefficient of 0.01 (and exploratory option of 0.03 to study the impact of higher coefficients); Use of the W-3 CHF correlation with non-uniform power distribution; Use of both Chen's [19] and Thom's [21] models of nucleate boiling, to compare the performance; Constant fuel gap conductance coefficient equal to 3070 W/m 2 K as estimated at average core burnup of 26.6 MWd/kgHM [1]; Temperature-dependent fuel and cladding thermal-physical properties [18].
Details of the sub-channel input model for VVER assemblies can be found in ref. [22].  Table 1 summarizes the suite of 13 variants of the MSLB scenario and the corresponding coarse-mesh computed core safety parameters. In each case, the parameters of the hottest assembly (which appears to be once burnt) and those of an adjacent assembly (twice or three times burnt) with high Fxy power peaking factor are shown. The nodal simulation results show that the highest return to power after scram and the highest radial peaking factors occur in Case 1 and Case 6. Figure 1 shows the time history of the total core power in Case 1 and Case 6 of the transient. The corresponding maximum return to power after scram is 1098 MW, respectively 1051 MW, at 67s from the onset of the transient. Figure 2 illustrates the computed relative assembly powers at time of max return to power (67s) in Case #1. The radial locations of the fully inserted control rod clusters are marked in blue and those of the stuck rods -in beige. The hottest assembly is the unrodded #104, marked in brick. Figure 3 illustrates the computed relative assembly powers at time of max return to power (67s) in Case #6. The core map format is as in Figure 2. The hottest assembly is the unrodded #129, marked in brick.

Full-core nodal calculations
The assembly-wise power distributions show that a big part of the core power is released in the region of stuck rods in the overcooled core sector. The peak of the 3D power distribution in the disturbed sector is near the core periphery. Correspondingly, Case 1, assembly #104 (15.3 MWd/kgHM) and the adjacent #117 (31.0 MWd/kgHM), as well as Case 6, assembly #129 (15.45 MWd/kgHM) and the adjacent #140 (24.1 MWd/kgHM) were selected for additional sub-channel thermal-hydraulic analysis. Figure 4 illustrates the coarse-mesh computed axial power distributions at time of max return to power -core-averaged and assembly-wise for the assemblies of interest for subchannel analysis.

Sub-channel thermal-hydraulic calculations
In Case 1, sub-channel TH analysis was carried out for the hottest assembly which is #104 (15.3 MWd/kgHM) and for the adjacent #117 (31.0 MWd/kgHM).
In both cases the assembly inlet and outlet TH boundary conditions at time of max return to power are nearly the same. The inlet temperatures for Case 1, assembly #104 and #117 are 207.62 ºC, and for Case 6, assembly #129 and #140 are 207.7 ºC and 207.66 ºC respectively. The corresponding inlet mass flow rates are nearly the same.
It should be noted that some preliminary sub-channel results for Case 1, hot assembly #104 have been reported in ref. [23], and this study extends the sub-channel analysis to more cases and various assemblies including such of higher burnup, where the UO2 melting temperature is lower.
The peak of the 3D power distribution in the disturbed sector is near the core periphery and corresponds to the radial location of the hottest assembly. For such a configuration and in the lack of full-core pin-by-pin calculation, it is reasonable to assume a hot full powerlike radial pin-power distribution for the hottest assembly (#104 in Case 1, and #129 in Case 6). An approximate tilt of ±5% across the assembly section is assumed for the adjacent assemblies. Figure 5 illustrates the adopted pin-power distribution at hot full power (HFP). Details of the computed distributions are illustrated in Figures 6 through 14. Table 2 summarizes the sub-channel results for the core safety parameters, obtained with the Chen and Thom options of the nucleate boiling model in CTF. The results in this table and the figures below (except for Figure 12) are obtained with user defined crosschannel mixing coefficient of 0.01. A small value is specified because this coefficient is not yet validated. Thus, the cross-channel mixing is mainly gradient driven. See Figures 11 and  12 for a parametric study of the impact of higher values of the user defined mixing coefficient on the vapor volume fraction.  In Table 2, the limiting values of the core safety parameters are as follows: Tmelt UO 2 , ºC: 2840 ºC for fresh fuel and 2540 ºC for burnt fuel (at 50 EFPD for the considered VVER fuel); Tclad: 1200 ºC; Min DNBR: 1.45 (when using the W-3 DNB correlation). The predicted values of the core safety parameters (within the limitations of the modeling assumptions) do not exceed the safety limits. However, the impact of the 5% tilt of the radial assembly pin power distribution is considerable and indicates that this issue requires further attention using full-core pin-by-pin coupled calculations. Figures 6 and 7 show the sub-channel results for Case 1, assembly #104 obtained with HFP pin-power distribution and Chen's nucleate boiling model. Figure 8 illustrates the results for Case 1, assembly #117, obtained with 5% tilt in the pin-power distribution across the hexagonal section and Chen's nucleate boiling model. Figures 9 through 12 show the results for Case 6, hot assembly #129 (15.45 MWd/kgHM). The results аrе obtained using HFP pin-power distribution. The comparison of the sub-channel liquid fraction maps in Figures 11 and 12 shows the impact of using higher values of the cross-channel mixing coefficient (0.03) -the sub-cooled boiling area shrinks toward the assembly center. Figures 13 and 14 show the results for Case 6, assembly #140 (24.2 MWd/kgHM), obtained with ±5% tilt in the pin-power distribution across the hexagon and Chen's nucleate boiling model.
The results in Figures 6, 8, 9 and 13 show that despite the high power peaking factors in the hottest assemblies only subcooled boiling and small vapor volume fraction in the upper part of the assemblies is predicted, because of the low inlet temperatures. Figure 9 shows that the use of Chen's nucleate boiling model predicts higher vapor void fraction in the subcooled region compared to that obtained with the Thom correlation. Both Chen's and Thom's correlations predict nearly the same axial profiles of the fuel centerline and cladding temperatures (see Figure 10).

Summary and conclusions
The safety parameters of a VVER-1000 core near the end of life were analyzed for a suite of pessimistic MSLB accident scenarios. For this purpose, coupled full-core COBAYA4/CTF neutronic/ thermal-hydraulic simulations at the nodal level were supplemented by CTF sub-channel calculations for the hottest assemblies of different burnup in the most risky cases.
In the CTF sub-channel assembly simulations, modeling options in the modeling assumptions with the Chen and Thom nucleate boiling models and with different userspecifies cross-channel mixing coefficients were compared.
The results show that: In all cases the predicted fuel cladding temperature at time of max return to power is far below the limiting value of 1200 ºC; The predicted sub-cooled boiling in the upper part of the hottest sub-channels is insignificant, even in the worst cases. The Chen model of nucleate boiling predicts higher vapor volume fraction compared to that obtained with Thom's model (as observed in other LWR assembly calculations); The minimum DNBR value of 2.66 predicted in the worst Case 1, Assembly #104 is well above the limiting value of 1.45 when using the W-3 CHF correlation; The predicted fuel centerline temperature in the hottest assemblies has a moderate margin of about 237 -251 ºC to the burnup dependent UO 2 melting temperature. Because of its importance, and of its strong dependence on the local pin power distribution and on the gap conductance coefficient, the assessment of this safety parameter deserves further attention.
In case of boiling in the assembly, the user-defined cross-channel mixing coefficient in sub-channel CTF needs further validation on measured assembly data because of its impact on the void fraction distribution. In the meantime, small values or not using this coefficient produce more conservative results. The results (within the limitations of the adopted modeling assumptions) show that the core safety parameters do not exceed the safety limits in the simulated aggravated reactivity accidents.
Further improved analyses are expected to include a coupled pin-by-pin neutronic/thermal-hydraulic calculation, at least for a mini-core of seven assemblies around the hottest one, as well as a dynamic gap conductance coefficient, a DNBR from VVERspecific CHF skeleton tables in CTF, and uncertainty quantification.