Abstract

In order to analyze the steady state and transient behavior of the CROCUS reactor, several methods and models need to be developed in the areas of reactor physics, thermal-hydraulics, and multiphysics coupling. The long-term objectives of this project are to work towards the development of a modern method for the safety analysis of research reactors and to update the Final Safety Analysis Report of the CROCUS reactor. A first part of the paper deals with generation of a core simulator nuclear data library for the CROCUS reactor using the Serpent 2 Monte Carlo code and also with reactor core modeling using the PARCS code. PARCS eigenvalue, radial power distribution, and control rod reactivity worth results were benchmarked against Serpent 2 full-core model results. Using the Serpent 2 model as reference, PARCS eigenvalue predictions were within 240 pcm, radial power was within 3% in the central region of the core, and control rod reactivity worth was within 2%. A second part reviews the current methodology used for the safety analysis of the CROCUS reactor and presents the envisioned approach for the multiphysics modeling of the reactor.

1. Introduction

A large variety of research reactors have been designed and operated during the last 50 years. These reactors are primarily designed for research purposes, yet they are widely applied in education and training, materials testing, and isotope production. Due to the diversity of research reactor designs and operating conditions, there is a wide variety of computational tools used in their safety analysis and, nowadays, it is desired to adopt a standard approach for safety analysis of these research reactors [1]. The development of high power research reactors and small modular reactors, together with the extended and intensive utilization of research reactors and the increased safety requirements of nuclear installations after the Fukushima accident [2], encourages the adoption of nuclear power plant (NPP) tools and methods to research reactor safety analysis. However, the use of NPP tools for research reactors is not straightforward as there are important differences in operating pressure, coolant flow, size, and power.

The coupling of thermal-hydraulic and neutronics codes becomes a fundamental tool for an accurate reactor behavior prediction under transient and accident conditions. Along those lines, a project financed by swissnuclear was started with the objective of developing methods and models for the coupled neutronics and thermal-hydraulics analysis of the CROCUS reactor at EFPL using advanced and state-of-the-art NPP computational tools. The present work represents the first stage of the project and focuses on the neutronics modeling.

This paper is divided into four sections. The first part briefly reviews the design of the CROCUS reactor. The second describes the methodology applied for the neutronics modeling of the reactor and the third section summarizes the process by which the model was benchmarked against a Monte Carlo solution. The fourth section reviews the current thermal-hydraulic modeling of the CROCUS reactor and describes the proposed model.

2. The CROCUS Reactor

The CROCUS reactor, operated by the École Polytechnique Fédérale de Lausanne (EPFL), Switzerland, is a two-zone uranium-fuelled, H2O-moderated critical research facility. It can be classified as a zero-power reactor, with a nominal power of 100 W. The core is approximately cylindrical in shape with a diameter of about 58 cm and a height of 100 cm. The reactivity in the CROCUS reactor is controlled by the water level, which can be adjusted with an accuracy of ±0.1 mm [3].

There are two different kinds of fuel rods within the CROCUS reactor core (see Figure 1). The central zone is fuelled with 336 UO2 fuel rods (1.806 wt%-enriched), which are thinner rods with a square lattice pitch of 1.8370 ± 0.0002 cm. The peripheral zone is loaded with 176 thicker, U-metal fuel rods (0.947 wt%-enriched) with a pitch of 2.9170 ± 0.0002 cm. All fuel rods have an aluminum cladding and are maintained in a vertical position by the upper grid and lower grid plates spaced 100 cm apart (see Figure 2). Because of the different pitches used, the two fuel zones are separated by a water gap, as shown in Figure 1. The core is located in an aluminum water tank of 130 cm diameter and 1.2 cm thickness. Light water (H2O) is used as moderator and reflector. With the current fuel loading, the critical water level is 95.22 ± 0.01 cm. Therefore, when the reactor becomes critical, a small axial section of the active core is exposed to air at atmospheric conditions as shown in Figure 2. There are two shutdown safety systems: expansion tanks that allow fast reduction of the water level and two cruciform control blades inserted from top to bottom. Figure 2 also provides a view of the reactor structure, the water tank, support plates, and fuel rods.

3. Neutronics Modeling

Although direct full-core transport calculations for transient analysis (such as DeCART [4], nTRACER [5], and MPACT [6]) are becoming possible with the increase of computational power, they remain very expensive and the full analysis of a nuclear reactor core currently relies on the traditional multistep methodology [7]. This approach begins with lattice physics to condense and homogenize spatially and spectrally the microscopic cross-section data into the structure needed for coarser-level codes (i.e., few-group parameters generation) and concludes with the core physics calculations to perform steady-state and transient full-core reactor calculations.

3.1. Cross-Section Generation

Traditionally, few-group parameters generation for full-core reactor simulators (such as PARCS) has been done using deterministic lattice physics codes. However, the use of continuous-energy Monte Carlo codes to generate few-group parameters can become an interesting option when dealing with reactor types that lie beyond the capabilities of conventional deterministic lattice physics codes [8]. CROCUS reactor characteristics make this methodology interesting as its core presents two incongruent fuel lattices with a water gap in-between, with no possible subdivision of the core in simple repeatable subsections (such as fuel assemblies).

Serpent 2.1.21, a Monte Carlo code developed at VTT [9], has been specifically designed for lattice physics applications. Serpent represents the state of the art for Monte Carlo lattice physics and has been chosen to provide the code PARCS with the homogenized cross sections. The use of Serpent code as a cross-section generator for PARCS code has been investigated by different research groups [10, 11].

In a previous work [12] the cross-section generation of the CROCUS reactor core was performed using Serpent code version 1.1.19. The SerpentXS python script [13] was used along with Serpent 1.1.19 to perform branch calculations and print cross sections into a PARCS compatible format. However, results from the previous work were not satisfactory since diffusion coefficients computed by Serpent 1 carried important errors of up to 30% [14]. The second release of the code, Serpent 2, implemented an updated approach for diffusion coefficient generation with improved accuracy [14].

Serpent 2 has the ability to generate diffusion coefficients using the classical definition based on the theory but it has also implemented fundamental mode methodology to correct diffusion coefficients based on an approximate leakage spectrum [8]. In this paper, both diffusion coefficients definitions were used and tested. The first definition, based on the traditional approximation, is computed aswhere is the microgroup diffusion coefficient, the macroscopic transport corrected cross section, the total macroscopic cross section, the average cosine of the scattering collision angle, and the zeroth moment of the scattering cross section. Then, the energy condensation (from micro- to coarse-group structure) of the diffusion coefficient is done as follows:where is the group index in the coarse-group structure.

On the other hand, when the leakage mode is invoked in Serpent 2, the code solves the equations [15] where is iterated to unity to get a better approximation of the neutron energy spectrum, resulting in the leakage-corrected flux spectrum and current spectrum . Then, the microgroup diffusion coefficient is computed aswhere is the energy independent buckling and is the group index in the microgroup structure. The energy condensation into a coarse-group structure is done with the leakage-corrected flux using (2).

Serpent 2 can solve the equations not only to provide an alternative definition of diffusion coefficient but also to generate leakage-corrected cross sections, that is, to use the critical flux spectrum for spectral collapsing of cross sections.

Since Serpent 2 uses different output variables names than Serpent 1 (the first release of the code), one of the tasks performed in this work was to update the SerpentXS scripts to become compatible with Serpent 2 output and to handle leakage corrected cross sections ( mode corrected). This updated script will be hereafter referred to as SerpentXS2.

Given that the CROCUS reactor core presents a peculiar geometry with two incongruent fuel lattices, the subdivision of the core in the form of fuel assemblies is not possible. For that reason, the most natural subdivision of the core is at a pin-cell level. This geometry was used to generate the fuel cross sections. Figure 3 illustrates the 2D geometry used to model the U-metal fuel (corresponding to the outer lattice) and also the one for the UO2 fuel (corresponding to the inner lattice). These heterogeneous pin-cell models use reflective boundary conditions in all three directions. Figure 3 also shows the difference in pin-cell sizes.

In order to generate the cross sections for control rods, a 2-D geometry of eight U-metal fuel rods with a control rod in the center was used as illustrated in Figure 4. Reflective boundary conditions were also used in all three directions. As shown in Figure 4, only the area surrounded by the dashed line was homogenized. The eight peripheral fuel pins were used to provide the heterogeneous problem with neutrons.

The water reflector region was modeled using a 2D geometry representing the radial boundary between core and reflector as illustrated in Figure 5. Reflective boundary conditions were used in the all directions with exception of the side of the reflector facing the vessel, which uses vacuum boundary conditions. Since in the bottom of the core there is a 47 cm layer of water, this heterogeneous model was also used to generate bottom reflector cross sections.

The top reflector region has been modeled with a geometry that included all structures on top of the core as shown in Figure 6.

Special treatment was taken over the water gap between fuel lattices to include it in either the UO2 or U-metal fuel lattices. For the case in which the water gap is contained in the outer lattice, the water gap volume was smeared across all U-metal nodes by increasing the fuel pitch from 2.917 to 3.023 cm. In a similar way, when the water gap is contained in the inner lattice, the UO2 pitch is increased from 1.837 to 1.909 cm.

For the cross-section generation of fuelled regions (i.e., UO2 and U-metal pin cells), two different diffusion coefficients definitions were used: a first one based on the approximation (1) and a second one based on the leakage model (3). Also, the cross-section spectral homogenization of all fuelled regions was computed in two different ways: one using the infinite flux spectrum (resulting from the infinite array of fuel pins) and a second one using the leakage-corrected flux spectrum (from the leakage correction model). Since the model is only applicable to regions where fission is taking place, reflector regions and control rods cross-section generation were limited only to -based diffusion coefficients and infinite flux spectrum for the cross-section collapsing.

Since the cross sections are generated at a pin-cell level, standard deterministic codes like CASMO [16], HELIOS [17], or TRITON [18] could have been used. Serpent 2 is chosen instead to take advantage of being able to model the full-scale heterogeneous problem, which represents the best available reference solution for the calculation scheme. A full-scale homogenization scheme is currently under development, which will provide an alternative to the presently used pin-cell level scheme.

Serpent 2.1.21 and the ENDF/B-VII nuclear data library were used for all Monte Carlo simulations (cross-section generation and full-core calculations). The full-core calculations were run using 900 cycles of 106 neutrons each, returning a final statistical uncertainty below 8 pcm for eigenvalue calculations and 0.1% for radial power distribution. For the pin-cell and other cross-section models, 1100 cycles of 105 neutrons each were used, returning a final statistical uncertainty below 0.01% for two-group parameters generation. The initial 100 cycles were skipped in all simulations.

3.2. Reactor Core Modeling

The first task of the reactor core modeling consisted of building a full-core model of the reactor using the Monte Carlo code Serpent 2. Since this model is used as reference for the comparison against the PARCS models, it was built including as many details as possible. In a previous work, the full-core Serpent model has been verified against a previously built MCNP model [19]. A validation work will be carried out in the near future.

Nodal methods are widely used for full-core reactor physics calculations. Each node normally corresponds to a small portion of the reactor core (e.g., to an axial slice of a fuel assembly) for which homogenized cross sections have first been obtained. PARCS [20] is a multigroup nodal diffusion code developed by the US NRC for 3D steady-state and transient analyses. However, PARCS also includes a finite difference kernel, which can be used for finer mesh solutions.

Due to the geometric characteristics of the CROCUS reactor core, subdivision of the core in fuel assemblies is not possible. Hence, two pin-by-pin full-core models were developed and run using the code PARCS v3.00. A first model uses a fine Cartesian mesh with a size equivalent to a UO2 fuel cell (1.837 cm) as illustrated in Figure 7. Since this finer mesh model cannot be used to predict power distribution in the outer lattice due to the mesh-fuel pin incongruences, a second model was required. The second PARCS model uses a coarser Cartesian mesh with a size equivalent to a U-metal fuel cell (2.197 cm) as illustrated in Figure 8. The latter can be used to predict the outer lattice radial power distribution; however, it fails to predict the inner lattice power distribution due to similar mesh-fuel pin incongruences.

PARCS calculations were run using two-group diffusion theory and a finite difference kernel. The two-group homogenized cross sections were generated using the Serpent 2 code as presented in the previous section. No correction factors such as interface discontinuity factors were used.

Axially, the active region of the core was subdivided into 25 nodes of 3.808 cm each, matching the 95.22 cm of water level as modeled in Serpent 2. Six nodes of 3.808 cm each were used to represent the region on top of the core.

Figure 9 shows the differences between the full-core Monte Carlo model and the two different PARCS nodalizations. The colors in the PARCS nodalization represent different cross-section sets: red for inner lattice, orange for outer lattice, green for control rods, and blue for reflector. Although it is difficult to visualize, the UO2 cell-size mesh model (Figure 9(b)) includes the water gap in the outer lattice cross-section set. Contrarily, the U-metal cell-size mesh model (Figure 9(c)) includes the water gap in the inner lattice cross-section set. The method used to include the water gap in one cross-section set or the other is described in Section 3.1.

4. Benchmark Results

This section is focused on the steady-state analysis and verification of PARCS results against a Serpent 2 full-core model. Three main steady-state parameters have been benchmarked: effective multiplication factor, control rods’ reactivity worth, and radial power distribution. The multiplication factor difference was computed as follows:where denotes each PARCS model.

The control rod reactivity worth was computed in Serpent 2 as the difference between a model containing the control rods fully withdrawn and the one with control rods fully inserted. In PARCS, the control rod worth was computed in a similar way, using a card that allows inserting or withdrawing the control rods. Finally, the percent difference reported in Table 1 was computed from the following expression:where denotes each PARCS model and stands for control rods.

The results shown in Table 1 suggest that PARCS models using -based diffusion coefficients and leakage-corrected cross sections are in good agreement with the Serpent 2 reference. The use of -based diffusion coefficients along with non-leakage-corrected cross sections results in eigenvalue underestimation. As for control rod reactivity worth, differences between the two PARCS models are related to the fact that control rods are being represented by three UO2 cells in the UO2 cell-size model and by only one U-metal cell in the other model (see Figures 9(b) and 9(c)). As a consequence, the UO2 cell-size model contains an 18% excess of control rod material.

Since a comparison does not return enough information on the overall accuracy of the PARCS model, an additional comparison exercise was performed focusing on radial power distribution. For this comparison, only -based models have been used. The two PARCS models, respectively, using the UO2 and U-metal cell-size mesh are used to predict the inner and outer lattice radial power distribution, since each one has a node-to-fuel pin matching in their respective regions as previously shown in Figures 7 and 8.

Figure 10 shows the steady-state radial power distribution predicted by Serpent 2 Monte Carlo code. The Serpent 2 radial power distribution was compared against the two PARCS models using the following expression:where the subindex denotes each node. Figure 11 illustrates the results of the radial power distribution comparison between Serpent 2 and PARCS. Note that the PARCS radial power distribution has been calculated using two different models: one for the inner lattice and another one for the outer lattice. The maximum nodal differences are of about 20% and are located in the inner lattice nodes adjacent to the water gap. Such differences could be related to the fact that both PARCS models smear the water gap, in the outer lattice for the UO2 cell-size mesh and in the inner lattice for U-metal cell-size mesh. PARCS power prediction in the central region of the core is within 3% with respect to Serpent 2. From safety analysis standpoint, these results are positive since the hottest rod will be most likely located in this area.

Since, as stated earlier, two different PARCS models were used to predict radial power distribution, additional verification was carried out between these two models and the Serpent 2 reference model. Table 2 shows the percent of the total power generated in each fuel lattice for the Serpent 2 model and also for the two PARCS models.

Table 2 differences between PARCS models and Serpent 2 are in the order of few percent. This could be potentially linked to the fact that while one PARCS model includes the water gap in the inner lattice, the other model includes it in the outer lattice.

5. Thermal-Hydraulic Modeling

For more than 50 years, numerous computer codes have been written to calculate the thermal-hydraulic characteristics of reactor cores under steady-state and operational transient conditions as well as hypothetical accidents. The main purposes of the continuing effort in the development of such computer codes have been improved computational effectiveness and improved ability to predict the response of the nuclear reactor. While thermal-hydraulic modeling plays a vital role in the design, operation, performance, and safety of a nuclear reactor, the present paper focuses particularly on the safety (or accident analysis) application.

Examples of codes used for research reactor thermal-hydraulic modeling are RELAP5 [21] for the NIST research reactor [22], the IPR-R1 TRIGA Brazilian reactor [23], and the High Flux reactor (HFR) in Netherlands [24]. The PARET code [25] has been used for the McMaster University research reactor [26], the University of Florida Training Reactor [27], and the NUR Algerian research reactor [28] among several others. Also, in many cases, thermal-hydraulic modeling was performed with in-house designed codes, such as the MULCH-II code for the MIT research reactor [29] and PLTEMP for the GRR1 Greek research reactor [30]. In all cases, the neutronic behavior of the research reactors was predicted using the point-kinetic approximation.

5.1. Current Methodology

The current safety analysis of the CROCUS reactor studies the reactor response under the maximum hypothetical accident, which is initiated by the flooding of the reactor core with light water at 12°C, the reactor being at nominal power (100 W), and failure of the shutdown safety systems. For the thermal-hydraulic modeling of this reactor, the in-house designed code EX_PUI [31] has been used. EX_PUI is based on a simple zero-dimensional model, assuming natural convection heat transfer and uniform temperatures in the fuel, moderator, and coolant. No two-phase flow equations are available in the code since the current accident analysis predicts that the maximum fuel cladding temperatures are of 60°C. The EX_PUI code also includes a point-reactor kinetic subroutine with reactivity feedback to predict the reactor power evolution. No overpower or power peaking factors are included in the current analysis.

The kinetic parameters used by the EX_PUI code were calculated using another in-house code, CRO93DIF [31], which is based on one-dimension multigroup diffusion theory. The CRO93DIF model of the CROCUS core used one-dimensional cylindrical coordinates. Four cross-section sets (for UO2 fuel, water gap, U-metal fuel, and reflector) and 19 energy groups were used for the diffusion calculations. The 19-group macroscopic cross sections used by CROF93DIF code were generated using the lattice code BOXER [32]. Also, the BOXER code was used to generate the reactivity feedback coefficients used by the point kinetic module of the EX_PUI code.

5.2. Proposed Methodology

The current CROCUS safety analysis relies on very simple models assuming a point reactor and all calculations are based on average values. A more detailed analysis, including multidimensional effects, power peaking factors, and the hottest channel analysis, may reveal restrictions or add flexibility to the reactor day-to-day operation. An additional driving force for this project results from the fact that the Swiss nuclear regulatory authority (ENSI) requested the Laboratory for Reactor Physics and System Behavior at EPFL, who is responsible for the CROCUS operation, to update the Final Safety Analysis Report using up-to-date tools and methods. Thus, the proposed update can provide additional details for the accident consequences and quantification of the conservatism in the original evaluations.

In general, many codes used for nuclear power plants can be also used for research reactor analysis. However, the ranges of parameters of interest to research reactors are different from those for nuclear power plants: this is namely true for fuel composition, system pressure, materials, and core geometry. Also, due to the large variety of research reactors, differences on validation and application may appear for each case. In this work, the proposed thermal-hydraulic analysis will be carried out using the TRACE code [33] since it is the current state-of-the-art tool for transient analysis of light water power reactors. It may be, however, necessary to modify TRACE to make it applicable for the particular channel geometry, coolant velocities, heat fluxes, and subcooled core conditions of the CROCUS reactor. As earlier described, research reactor thermal-hydraulics models have been coupled to point-reactor kinetic modules [34]; however, the proposed method takes advantage of the relatively simple coupling of TRACE to PARCS for the multiphysics modeling. By doing so, we are not only utilizing state-of-the-art methods for research reactor analysis but also investigating TRACE/PARCS potential applications to small modular reactors (SMR).

In a first approach, the TRACE model of the CROCUS reactor will consist of a 3D vessel component for the reactor vessel. The core will be represented by several heat structures that will be coupled to the PARCS model, which solves the neutron kinetic problem. Thus, power evolution, peaking factors, and reactivity feedback coefficients are computed by PARCS and transferred to TRACE. Several hydrodynamic channels will be used within the vessel component to represent the different areas of the core. Reactor parameters and operating conditions considered in the safety analysis will be chosen assuming the most unfavorable conditions, that is, following a conservative approach.

6. Conclusions

In this paper, a methodology for the coupled neutronics and thermal-hydraulics analysis of the CROCUS reactor at EFPL has been studied. The Serpent 2 Monte Carlo code has been used to generate two-group parameters for the PARCS code. Since the Monte Carlo technique offers significant advantages for detailed modeling of the complex geometrical configuration of the reactor core, a full Serpent 2 model of CROCUS has been built. A steady-state benchmark has been conducted between PARCS and Serpent 2 full-core models. Good agreement was achieved in terms of eigenvalue calculations and control rod reactivity worth. Radial power distribution results show good agreement in the central core region; however, they reveal that PARCS models present limitations in predicting power near the water gap region.

Future work will address validations of the PARCS model for the static and dynamic analysis. Also, since the full-scale heterogeneous problem represents the best available reference solution for cross-section homogenization, it is currently under development to provide an alternative and a potential improvement over the pin-cell homogenization scheme. Future tests may also include the use of superhomogenization factors (SPH) [35] to yield better approximation of the full heterogeneous problem.

The PARCS model will be coupled to a TRACE thermal-hydraulic model of CROCUS for transient analysis of the reactor. From an accident analysis perspective, the current CROCUS safety analysis report shows room for improvement as it relies on very simple models that may unnecessarily limit the range of operation of the reactor. Thus, reassessment using state-of-the-art tools would provide not only more realistic predictions that reduce the deliberate conservatism but also the possibility to add flexibility in the day-to-day reactor operation.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.