Next Article in Journal
Model Predictive Control of a Modular 7-Level Converter Based on SiC-MOSFET Devices—An Experimental Assessment
Next Article in Special Issue
A Qualitative Strategy for Fusion of Physics into Empirical Models for Process Anomaly Detection
Previous Article in Journal
Multi-Objective Constructal Design for Square Heat-Generation Body with “Arrow-Shaped” High-Thermal-Conductivity Channel
Previous Article in Special Issue
Projecting the Thermal Response in a HTGR-Type System during Conduction Cooldown Using Graph-Laplacian Based Machine Learning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

CTF-PARCS Core Multi-Physics Computational Framework for Efficient LWR Steady-State, Depletion and Transient Uncertainty Quantification

Department of Nuclear Engineering, North Carolina State University, Raleigh, NC 27695-7909, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Energies 2022, 15(14), 5226; https://doi.org/10.3390/en15145226
Submission received: 1 April 2022 / Revised: 26 April 2022 / Accepted: 14 July 2022 / Published: 19 July 2022
(This article belongs to the Special Issue Latest Advances in Nuclear Energy Systems)

Abstract

:
Best Estimate Plus Uncertainty (BEPU) approaches for nuclear reactor applications have been extensively developed in recent years. The challenge for BEPU approaches is to achieve multi-physics modeling with an acceptable computational cost while preserving a reasonable fidelity of the physics modeled. In this work, we present the core multi-physics computational framework developed for the efficient computation of uncertainties in Light Water Reactor (LWR) simulations. The subchannel thermal-hydraulic code CTF and the nodal expansion neutronic code PARCS are coupled for the multi-physics modeling (CTF-PARCS). The computational framework is discussed in detail from the Polaris lattice calculations up to the CTF-PARCS coupling approaches. Sampler is used to perturb the multi-group microscopic cross-sections, fission yields and manufacturing parameters, while Dakota is used to sample the CTF input parameters and the boundary conditions. Python scripts were developed to automatize and modularize both pre- and post-processing. The current state of the framework allows the consistent perturbation of inputs across neutronics and thermal-hydraulics modeling. Improvements to the standard thermal-hydraulics modeling for such coupling approaches have been implemented in CTF to allow the usage of 3D burnup distribution, calculation of the radial power and the burnup profile, and the usage of Santamarina effective Doppler temperature. The uncertainty quantification approach allows the treatment of both scalar and functional quantities and can estimate correlation between the multi-physics outputs of interest and up to the originally perturbed microscopic cross-sections and yields. The computational framework is applied to three exercises of the LWR Uncertainty Analysis in Modeling Phase III benchmark. The exercises cover steady-state, depletion and transient calculations. The results show that the maximum fuel centerline temperature across all exercises is 2474 K with 1.7 % uncertainty and that the most correlated inputs are the 238U inelastic and elastic cross-sections above 1 MeV.

1. Introduction

Best Estimate Plus Uncertainty (BEPU) [1] approaches have received a lot of interest during the last few years in the context of nuclear reactor safety analysis. Efforts are being made to improve the multi-physics and multi-scale computational modeling of nuclear reactors. Various coupling approaches have been used to establish robust computational models under different conditions from steady state to transient with reasonable accuracy and computational efficiency. Typically, such modeling approaches restrain the multi-physics analysis to a coupling between thermal-hydraulics and neutronics. The fuel performance is simplified to capture the main feedbacks through the thermal solver inside the thermal-hydraulics modeling. This is performed because fuel performance is computationally very expensive due to the complexity of the multi-scale phenomena involved. Various computational frameworks have been developed and established around the world to develop such multi-physics calculations. In the U.S., VERA [2], NEAMS [3] and MOOSE [4] are the main examples that can have an even higher multi-physics fidelity. In France, CEA has developed its own platform using SALOME [5], while EPFL and PSI in Switzerland have explored the OpenFOAM capabilities to establish their platform [6]. In Finland, VTT has recently developed a new multi-physics platform called Kraken [7] and in Korea, UNIST has developed MPCORE [8], a multi-physics coupling that can even model the core at a pin-by-pin scale. In this work, we present the core multi-physics computational framework for efficient Light Water Reactor (LWR) steady state, depletion and transient uncertainty quantification established at North Carolina State University (NCSU) [9]. In this framework, the subchannel thermal-hydraulic code CTF [10] and the nodal expansion neutronic code PARCS [11] are coupled at the reactor core level. This CTF-PARCS coupling leverages the ongoing developments of CTF at NCSU.
A wide variety of methods are available for coupling the different physics. A first classification of these coupling methods is usually made by dividing them into tight and loose coupling [12]. In tight coupling, the different physics are embedded in a unique code that solves all the constituent equations. On the other hand, in the loose coupling, the different physics are solved separately using an external iteration scheme to track the convergence of the coupled solution. Furthermore, loose coupling methods can be divided into internal and external coupling schemes [13]. The internal coupled codes are compiled in the same project, and the communication and data exchanges take place by sharing their internal memory. On the other hand, the externally coupled codes remain isolated in their own compilation projects, and the communication and synchronization relies on a coupling interface based on any of the available inter-process communication protocols. Different solution strategies are available for solving the coupling of the different physics and their effectiveness depends on the above-mentioned coupling methods. For example, for tightly coupled codes, the implementation of the implicit Jacobian-free Newton–Krylov (JFNK) method [12] is desirable because of its high convergence rate. Intrusive and extended code modifications are needed for the implementation of the JFNK method. For this reason, although conceptually it can be applied in loosely coupled codes, it is usually reserved for tight coupling. For loose coupling methods, the Picard iteration method [14] is widely used due to its simplicity. This method has the associated problem of poor convergence rate, but this can be remedied by applying acceleration methods such as the Alternating Nonlinear method, the Residual Balance, or the more convenient Adaptive Residual Balance [15]. The Semi-Implicit Operator Splitting [12] is another well-developed and optimized solution strategy to consider in loosely coupled codes. Applications of these methods and solutions strategies to the neutron kinetics and thermal-hydrualics coupling can be found in these references [16,17,18,19]. From the presented methods and solution strategies the loose coupling and the Picard iteration have been selected for the CTF-PARCS coupling. No acceleration methods are currently implemented apart from the relaxation of the relative tolerance of CTF.
The developed CTF-PARCS coupling is embedded in an uncertainty quantification framework that allows efficient and robust statistical analysis of various input and output quantities of interest. This allows the framework to be applicable to BEPU approaches that alleviate some of the conservatism involved in safety studies. In this context, the OECD/NEA Light Water Reactor Uncertainty Analysis in Modeling (LWR-UAM) Benchmark [20] was initiated to facilitate the development of uncertainty quantification methodologies and propagate consistently uncertainties across multi-scale and multi-physics coupled calculations. The application of BEPU approaches to transient core calculations is challenging due to the increased computational cost. In [21], the CEA multi-physics platform was used to perform such studies, while in [22], a coupling between PARCS and system thermal-hydraulics code TRACE [23] was used for I-2a, I-2b and I-2c Pressurized Water Reactor (PWR) exercises of the LWR-UAM Phase III [24]. In this work, we perform a similar uncertainty quantification study using the developed CTF-PARCS computational framework for the same exercises of the LWR-UAM Phase III. These exercises involve the modeling of the TMI-1 core under various conditions. For Exercise I-2a, we model the steady-state Hot Full Power (HFP) at the Beginning of Cycle (BOC) and End of Cycle (EOC). For Exercise I-2b, we perform the depletion as a multi-state coupled calculation from BOC to EOC. For Exercise I-2c, we model a rod ejection accident (REA) initiated from HFP conditions at both BOC and EOC.
This article is structured in five sections; Section 1 was the above presented introduction. Section 2 discusses in detail the developed multi-physics computational framework. The different codes that are being coupled are first presented followed by a description of the Python tool used to streamline the transition from lattice to core neutronics. The implemented CTF-PARCS coupling approaches for steady-state, depletion and transient calculations are discussed, and the uncertainty quantification approach is presented. Section 3 describes the LWR-UAM Phase III uncertainty quantification exercises selected to apply the developed computational framework. Section 4 presents and analyzes the obtained results for all the studied exercises. Finally, in Section 5, the general conclusions from this work and future intended developments are discussed.

2. Computational Framework

This section describes the computational framework used to simulate the proposed steady-state, depletion and transient calculations. It focuses on the capabilities and roles of the standalone codes involved in the framework, as well as the interaction of the codes embedded in the framework to perform the calculations and the adopted uncertainty quantification approach.

2.1. Simulation Tools

This subsection contains a short description about each one of the simulation tools involved in the developed computational framework. This includes the main capabilities and the most important features of these codes used in our current activities.

2.1.1. Polaris

Polaris [25] is a two-dimensional (2D) deterministic neutron transport code from the SCALE 6.2.4 suite [26]. Polaris is used to produce multi-history, multi-branch cross-section libraries at different burnup. The cross-sections are condensed from one of SCALE 6.2.4’s predefined fine-group structure into a coarse-group structure. The generation of the cross-sections follows a traditional approach, where each fuel assembly representative of the LWR core of interest is modeled individually in Polaris to generate a set of homogenized few-group cross-sections in the SCALE format. Reflective boundary conditions and a critical spectrum are applied in these lattice models. Polaris can also generate homogenized few-group cross-section libraries for the reflector regions of the core. The models implemented to generate the nominal set of cross-sections are described in Section 3.2.
The generation of the two-group cross-sections with Polaris results in a set of SCALE-formatted libraries (hereafter referred to as t 16 format). Each library contains cross-sections representative of one fuel or reflector assembly constitutive of the core. The Polaris libraries contain information at different burnup and branch points. The burnup points are required in the cross-section modeling to account for the burning, breeding, decaying and activation of actinides and fission products found in the materials depleted in the Polaris models. Branches are used to model instantaneous state variations (temperature, density, soluble boron concentration, etc.) during operations and transients. Histories account for the long-term effects of the neutron spectrum on the burned materials. The historical parameter is modeled separately in Polaris, which means that modeling H different historical points for a given fuel assembly requires H Polaris inputs. Overall, a core containing A F types of fuel assemblies and A R types of reflector assemblies with H F fuel assembly histories and H R reflector assembly histories requires T = A F × H F + A R × H R Polaris inputs and produces T libraries in t 16 format.

2.1.2. Sampler

Sampler is a module of SCALE 6.2.4 [26] used in association with Polaris to generate perturbed cross-section libraries. Sampler is used in the framework to perturb the multi-group cross-sections, the fission yields of fissionable actinides, and manufacturing parameters of the lattice cell. Perturbing consists of modifying simultaneously the uncertain inputs (cross-section, yields, manufacturing parameters) in the assembly models to obtain sets of cross-sections libraries representative of the input perturbations. The perturbed cross-section libraries are used in the CTF-PARCS simulations to assess the propagation of uncertainties (cross-sections, yields, manufacturing parameters) into key core outputs of interest.
In Polaris-Sampler, the range and nature of the sampling distribution for the manufacturing perturbations are user-defined and are described for our application in Section 3.2. The number of degrees of freedom is limited to the number of samplings for the cross-sections and fission yields perturbations, because perturbation factors have been pre-computed in the SCALE libraries. The perturbation factors were generated in SCALE from the covariance data available in ENDF/B-VII.1 [27] along with low-fidelity covariance data supplemented for minor fission products and actinides.

2.1.3. GenPMAXS

The Polaris generated libraries in t 16 format are not readable by the core simulator, PARCS (Section 2.1.4), and need to be converted into a PMAXS format using the GenPMAXS sequence [28].
The GenPMAXS program combines into a PMAXS library the history, branch and burnup points of a given assembly. The construction of the GenPAMXS inputs relies on:
  • The configuration of the fuel or reflector assembly defined by the assembly specifications, such as the pin pitch, number of pin per assembly. Those specifications are required in the PMAXS for pin power reconstructions.
  • The cross-section parameterization (branch, history, burnup). GenPMAXS also requires the definition of a reference history and branch.
  • Pre-defined user inputs. For example, the name given to the PMAXS file, boolean flags such as inclusion of assembly and axial discontinuity factors.

2.1.4. PARCS

PARCS is a 3D reactor core neutron kinetics simulator code that solves the steady-state and time-dependent multi-group neutron diffusion and low-order transport equations (SP3) in Cartesian and hexagonal geometries [11]. Two additional important features of PARCS are the capability of reading burnup-dependent cross-sections in PMAXS format and making pin power reconstruction when group-dependent flux form functions are provided within the cross-sections library. PARCS exists as a standalone code and as a coupled version to TRACE. The TRACE-coupled internal version is bundled directly into the TRACE distribution and has been used in this project as it facilitates the multi-physics multi-scale coupling between TRACE, CTF, and PARCS.

2.1.5. CTF

CTF is the shortened name given to the version of COBRA-TF being developed and improved by the Reactor Dynamics and Fuel Modeling Group (RDFMG) at North Carolina State University (NCSU) [10]. CTF is a subchannel thermal-hydraulic simulation code designed for Light Water Reactor (LWR) vessel and core analysis. It uses a two-fluid three-field modeling approach for the two-phase flow. In the last decade, CTF has been extensively developed and validated for PWR, BWR, VVER, Small Modular Reactor (SMR), Fast Breeder Reactor (FBR), and research reactor applications. CTF is a state-of-the-art subchannel code for reactor thermal-hydraulics bundle and core analysis and is a part of the U.S. DOE CASL [29] and E.C. NURESAFE [30] projects.

2.1.6. Dakota

Dakota [31] is a general toolkit to perform uncertainty quantification, optimization and sensitivity analysis for any code that can be treated as a black box. In this work, the code is the multi-physics coupling between CTF and PARCS, and Dakota is used for performing the uncertainty propagation and more specifically the sampling of CTF input parameters and boundary conditions and the parallel distribution and execution of the different calculations. Custom Python scripts are used for providing an interface between the code and Dakota and for post-processing and analyzing the obtained results.

2.2. Polaris Pre-Processing

In general, 50 to 100 Polaris inputs are required to construct a full-core LWR model in PARCS to account for the core history of each assembly type in nominal conditions. The construction of multi-history inputs of different assemblies requires the modification of the following parameters in the Polaris inputs:
  • The pin layout.
  • The material definitions.
  • The reference state.
  • The core plate thickness or shroud thickness, if any, for reflector assemblies.
The historical, branch and burnup points can remain identical for each fuel assembly. The same statement applies for reflector assemblies, although the reflector parameterization does not have to align with the fuel assemblies’ parameterization. After establishing the parameterization of the nominal configuration of each assembly, Sampler creates N perturbed cases per assembly history. The execution of Polaris/Sampler produces N t 16 libraries. These libraries must be converted into the PMAXS format; therefore, the GenPMAXS inputs must be implemented for each assembly required by PARCS for each perturbation computed used for the uncertainties analysis. The manual implementation of the Polaris and GenPMAXS inputs is error-prone, redundant and time-consuming; therefore, the generation of these inputs as well as the code executions were streamlined into a unique data management sequence. The sequence accelerates the generation of the PMAXS libraries, which enhances the consistency and robustness of the libraries. It relies on two Polaris templates files: one for the fuel assemblies and one for the reflector assemblies. Both templates are formatted in the Sampler formulation to avail the generation of the perturbed assemblies. One of the modules of the sequence contains the properties of the assemblies, such as the pin layouts of fuel assemblies and the historical points chosen for the cross-section parameterization. The branch points and the burnup points are hard-coded into the templates because they are uniform across all the fuel and reflector assembly models. Optionally, the dimension and composition of the pins can be implemented in the sequence with corresponding uncertainties, if manufacturing uncertainties are included in the model. The construction of the Polaris assemblies is illustrated in Figure 1.
After execution of the sequence, the creation of the Polaris inputs is separated into two phases: (a) the nominal values are used to produce the nominal Polaris assemblies, and (b) the Sampler inputs are built; therefore, the Sampler aliases are preserved from the template, and a distribution block is added at the end of the Polaris inputs using the nominal values and the associated standard deviations. After the construction of the nominal and perturbed inputs, the nominal cases are executed, which are followed by the perturbed cases. Typically, 100 libraries per assembly and per history leads to about 10,000 perturbed Polaris inputs/outputs. If the computational power available is sufficient to run the entire set of nominal and perturbed Polaris simulations, the sequence will submit the jobs and terminate when the t 16 libraries exist. Otherwise, the sequence can be stopped while the Polaris simulations are running and be resumed at a later time. The sequence is set to submit the simulations on a Linux High-Performance Computing with a Slurm workload; the portability of the sequence depends on the Operating System and the workload. Additionally, the perturbation of manufacturing parameters involves the perturbation of the pin dimensions or reflector dimensions in the assembly, which can induce ray spacing errors in SCALE 6.2.4. It was found that ~8% of the Polaris inputs failed due to ray spacing error with the default spatial discretization; this number tends to increase if the spatial discretization is refined. To handle robustly the situations of (a) resuming the submission of the simulation on the HPC or (b) fixing ray spacing issues, the sequence reads the existing Polaris outputs to re-run or re-index cases. Therefore, a set of M perturbed inputs are generated, while N perturbed are actually run, to have E inputs cases ready to substitute the failed input configurations (see Equation (1)).
E = M N , M > N , { N , M } N
At a given perturbation n , 1 n N , the sequence restarts are handled as follows:
  • If a SCALE execution file exists ( e x e c file), the job is currently running, and the sequence moves on to the next perturbation, n n + 1 .
  • If no t 16 file and no Polaris output ( . o u t ) exist, the perturbation case has not been run yet: the sequence runs the simulation normally and moves on to the next perturbation, n n + 1 .
  • If the t 16 file and Polaris output file both exist, the number of state points computed in the t 16 file ( S t 16 ) is cross-compared to the expected number of state points predicted by the Polaris output ( S out ).
    -
    If S t 16 = S out , the perturbation case has already run normally, and the sequence moves on to the next perturbations, n n + 1 .
    -
    Otherwise, the sequence raises an error, and no diagnostic is predicted.
  • If the Polaris output exists, but the t 16 does not, the Polaris output is read by the sequence to find out if the simulation stopped due to a ray spacing error.
    -
    If a ray spacing error occurred, another perturbation is drawn from the E supplemental inputs. The supplemental case indexed m is re-indexed into n , the case is re-run, and then the sequence moves on to the next perturbations, n n + 1 .
    -
    Otherwise, a job crash is assumed (exceeded wall time, reboot, etc.) and the perturbation n is re-run without modification and then moves onto the next perturbation, n n + 1 .
A summary of the execution and verification phase is provided in Figure 2. Once all the nominal and perturbed t 16 libraries are generated, the sequence advances to the PMAXS generation phase.
The GenPMAXS files are constructed immediately from the information contained in the Polaris input files without using templates. The type of reactor, the name of the assembly, the number of history points, the number of energy groups, the types of state parameters (boron concentration, fuel temperature etc.), the values of the historical points, the number of historical points, the number of branches, the values of the branch points, the number of burnup points, the values of the burnup points and the names of the t 16 libraries are extracted from the Polaris input. The sequence selects a reference history point based on the type of reactor by calculating the Euclidian distance between each historical points to the predefined reference conditions of the reactor. The reactor conditions involve the fuel temperature ( T F ), the moderator temperature ( T M ), the moderator density ( ρ M ) and the boron concentration ( C B ). For example, the PWR reference conditions are defined in the sequence with T F = 900 K, T M = 573 K, ρ M = 0.726 g · cm 3 , C B = 900 ppm. This approach is repeated to determine the reference branch point. The same approach also operates for the reflector, with the exception that PMAXS libraries are generated for corner reflectors in addition to the regular reflector libraries. Corner reflector libraries are computed from regular reflector libraries, using an additional parameter r 2 D in GenPMAXS computed from the assembly pitch P FA and the shroud thickness d , as shown in Equation (2).
r 2 D = P FA d P FA
After the generation of the GenPMAXS inputs, the sequence executes automatically the GenPMAXS inputs to generate the PMAXS data for the both the nominal case and the perturbed case.

2.3. CTF-PARCS Coupling

The subchannel thermal-hydraulic code CTF and the core neutron kinetics simulator PARCS have been recently coupled at NCSU. The capabilities of this coupled code are being expanded, and the code can be currently used as a three-dimensional multi-physics core simulator for steady state, depletion and transient. This core simulator, referred to as CTF-PARCS in this manuscript, will be used in the multi-physics simulation framework for uncertainty quantification proposed here. A flexible Message Passing Interface (MPI) communication protocol based on a server–client coupling algorithm has been developed as a kernel of our multi-physics multi-scale simulation platform. CTF-PARCS is the core simulator in the general simulation platform developed at NCSU and presented in Figure 3. This coupling algorithm facilitates the implementation of loose coupling methods between different simulation tools without the necessity of extended modification in their sources and compilation projects [32].
In the developed CTF-PARCS multi-physics core simulator, CTF acts as the server while PARCS is the client, as indicated in Figure 4. The coupled simulation is activated by using a command line argument in CTF, being as a server able to launch additional MPI processes acting as clients, as will be PARCS in this case. A coupling interface input file (input-name.ci) handles the information about the different clients launched by the server. The neutron kinetics coupling class developed in CTF contains protocols to control the run-time advance of PARCS, such as initialization, steady-state iteration, depletion step, transient time step, convergence, finalization, and error. The synchronization is managed through tags communicated using the MPI standard functions MPI_SEND and MPI_RECV. The feedback parameters (included in Figure 5a) are exchanged during each global coupled code iteration using the MPI communication PORT opened between both codes. Mapping and auto mapping procedures have been developed to provide the information about the different physical domains, spatial discretization, and spatial mesh overlays. It must be noted that CTF is a MPI parallel code, and different approaches can be taken regarding the subdomain decomposition. Currently, the MPI process with rank 0 in CTF is making the full domain communication in the coupling algorithm following in fact a global mapping scheme provided as input to CTF.
To analyze the global convergence of the core simulator, the three developed simulation modes, steady state, transient, and depletion, must be analyzed separately.
  • In steady-state mode, a simple Picard iteration scheme [14] has been adopted. This method implies the construction of an outer convergence loop to check for global convergence of the coupled code and as many local convergence loops as physics are being coupled together. These nested loops represent the global coupled convergence. In this method, the convergence of the thermal-hydraulic code has to be fully met before providing feedback to the neutronic code. The neutronic code then uses the thermal-hydraulic feedback to resolve the neutron flux distribution until its own convergence criteria are met. The global coupled code convergence is tracked by the residuals of the feedback parameters using the L 2 and L norms of each parameter until the defined tolerance is reached. The L 2 and L norms are defined in Equations (3) and (4), respectively, where i is the cell number, nc is the total number of cells, j is the index of the residual of a given feedback parameter, and nr is the total number of residuals. R i , j is the residual of the j th feedback parameter in the i th cell for the current iteration, which is computed as R i , j = R i , j n R i , j n 1 , among which n is the current iteration number. The coupled code steady-state convergence behavior is illustrated in Figure 6.
    R 2 = i = 1 nc j = 1 nr | R i , j | 2
    R inf = max ( | R i , j | i = 1 , , nc j = 1 , , nr )
  • For transient simulations, a temporally explicit coupling scheme has been adopted. Therefore, the neutronic and thermal-hydraulic conservation equations are not solved altogether, but each code solves its own system of transient equations using the most limiting time step size of both codes. This explicit client/server coupling scheme is presented in Figure 5b. The temporally explicit methods are usually adopted when coupling thermal-hydraulics with neutron kinetics because the numerical stability is always preserved. Although it is widely adopted, the temporally explicit coupled methods have two main inconveniences. The first issue is the numerical diffusion, as each of the coupled physics relies on the solution of the other physics at the previous time step. The second inconvenience is the limitation of the time step size, as all the coupled codes must adopt the smaller time step size between all the involved physics.
  • In depletion mode, a multi-state simulation procedure is initialized in the core simulator. The number of states of this multi-state mode is defined by the number of depletion steps set in PARCS input. At the beginning of each state, the 3D burnup distribution is sent from PARCS to CTF to initialize the burnup-dependent models and material properties in the fuel rod objects. After this initialization, the Picard iteration global convergence scheme is used in the same way as in steady-state simulations. It must be noted that the boron concentration is not feedbacked from CTF to PARCS in depletion mode, as the critical boron is the parameter varied by the latter to find the criticality of the core. This multi-state iterative procedure is represented in the scheme of Figure 7. The global convergence of the developed core simulator in depletion mode can be analyzed by looking at the L norms presented in Figure 8. In Figure 8, it can be observed that the global convergence is longer in the first iteration because the initial conditions are far from the real solution. All the remaining states converge with a similar convergence ratio. Each peak on the residual represents the beginning of each state. The peaks are large due to the change in the 3D burnup distribution in the core.
Some additional capabilities have been implemented in CTF during the development of CTF-PARCS coupling for improving the prediction of the effective fuel temperature:
  • The first development is the possibility of defining customized ways to compute the Doppler fuel temperature. This effective fuel temperature can be defined as a linear combination of the pellet radial average temperature ( T avg ), the fuel surface temperature ( T surf ), and the fuel centerline temperature ( T cen ), as shown in Equation (5). The parameters A , B , and C are defined by the user in the coupling interface input file.
    T eff = AT avg + BT surf + CT cen
    This allows the user complete flexibility on the definition of the effective temperature. Typical effective fuel temperature formulas found in the literature that can be reproduced with this development are the Rowlands [33], as shown in Equation (6), and the Santamarina [34], as shown in Equation (7).
    T eff = 4 9 T cen + 5 9 T surf
    T eff = T avg 1 18 ( T cen T surf )
  • The second added feature is the initialization of the pellet power and exposure radial profiles in CTF. This capability has been implemented by linking an external 1D depletion module to CTF. This external depletion module can be called from the fuel rod class of CTF to make a depletion and obtain the radial power and exposure profiles for each nuclear rod and axial node of the model.
In this project, the effective fuel temperature is calculated using the Santamarina formula [34] of Equation (7), and the radial power and exposure profiles within fuel pellets are calculated using the above-mentioned capability. For the latter, we have selected the TUBRNP model developed in TRANSURANUS code [35] from the publicly available model to compute the radial dependence of the power and burnup in the fuel pellets. This model determines the power profile based on the fuel rod geometry and the initial fissile material concentration. Single group thermal flux diffusion theory is applied to calculate the concentration of 235U, 238U, 239Pu, 240Pu, 241Pu, and 242Pu isotopes in each pellet ring. The microscopic cross-sections are spatially constant. Bessel functions are used to solve the 1D single group diffusion equation, obtaining the flux shape within the pellet as a function of the isotopic composition of each ring. For further details, the reference [35] can be consulted. This model has been coded and compiled as an external library and is used to initialize the fuel rod power and burnup profile at the beginning of the CTF-PARCS steady state and before every one of the depletion steps. PARCS sends the 3D burnup distribution to CTF at the beginning of the simulation and at the end of each depletion step.

2.4. Uncertainty Quantification

The previously discussed multi-physics computational model is embedded into an uncertainty quantification framework, as illustrated in Figure 9. The framework combines different custom Python scripts with Dakota functionalities to perform a consistent and efficient uncertainty quantification study. The consistency of the framework requires applying all the multi-physics inputs perturbations to all the relevant physical domains. For example, if the fuel rod radius is perturbed, this perturbation needs to be applied to Polaris but also to be considered in the fuel rod geometry, in the calculations of the flow areas, and in the wetted perimeters in CTF. The custom-made Python scripts allow an easy and straightforward application of various multi-physics inputs. Once the perturbations are applied, Dakota distributes the calculations in parallel and executes them until all the outputs have been obtained. The sampled inputs and outputs are used for the statistical analysis, where different statistical quantities such as mean, standard deviations and correlations are calculated.
The work emphasizes estimating the correlations between multi-physics outputs and the 56-group microscopic cross-sections available in SCALE and used by Sampler and Polaris. In state-of-the-art applications, the two-group macroscopic cross-sections are one of the figures of merit correlated to the downstream outputs. The objective here is to (a) retrieve a finer level of information than two-group macroscopic data, and (b) correlate the uncertain inputs from a more upstream level (i.e., the 56-group cross-sections before the lattice calculations) to the multi-physics outputs. The fundamental difficulty in such conventional two-step multi-physics calculations is that the computational cost does not allow any detailed sensitivity analysis to determine the important inputs. Recent studies [22] addressed this issue by perturbing independent separate groups of inputs and evaluating the standard deviations obtained for each group with the standard deviation obtained when all the inputs are perturbed altogether. The issue with such approaches is that (a) it requires larger (but feasible) number of calculations, and (b) the grouping of inputs is done at a very coarse level. For example, the inputs are separated by physics domains (e.g., neutronics, thermal-hydraulics). The large number of inputs that are highly correlated constitutes the major challenge for going to a finer level of sensitivity analysis. This is specifically true for neutronics due to the large number of microscopic cross-sections and yields stored in the continuous energy and multi-group libraries. The other sources of input uncertainties (e.g., manufacturing parameters, boundary conditions) can be assumed independent with reasonable assumptions and do not raise problems for the statistical analysis. In the developed computational framework, we exploit the existence in SCALE of pre-calculated perturbation factors for all the 56-group microscopic cross-sections and yields. The PALEALE utility of the SCALE suite was used to pre-process and extract all these perturbation factors. This pre-processing needs to be done only once, since the perturbations are constant and they apply to every uncertainty quantification study performed with Sampler. The obtained microscopic cross-sections and yields perturbations can then be used to estimate their correlations with the multi-physics outputs of interest. Due to the large quantity of data, only a selection of the most important isotopes is analyzed: 235U, 238U, 239Pu, 1H, and 16O. For these isotopes, the reactions listed in Table 1 are considered.
Although the yields have been studied, their impact in this work is found negligible for the analyzed cases and thus is not discussed. The Spearman rank correlation coefficients are used for estimating the correlations. Spearman coefficients are computed based on the ranking of the variables and thus capture monotonic trends. If we assume a sample of size N where R x and R y are the vectors of ranks for each sample of two quantities x and y , respectively, then the Spearman coefficient r xy can be computed through Equation (8). The Spearman coefficients are preferred over the typical Pearson correlation that captures only linear relationships. In the future, more advanced methods such as the Hilbert–Schmidt Independence Criterion (HSIC) indices could be considered [36].
r xy = 1 6 i = 1 N ( R x i R y i ) 2 N ( N 2 1 )
The developed statistical framework allows the calculations of correlations between multi-physics outputs and the 56-group SCALE microscopic cross-sections. Before presenting the different studies and their results, it is important to discuss three aspects relevant to such an approach:
  • The first aspect is the statistical noise. Due to the limited number of samples, even independent quantities can be correlated up to some degree. For this reason, a cutoff threshold needs to be determined in order to screen out all the correlations that are within the statistical noise. To calculate the cutoff threshold, a very simple brute force approach is performed. For a predefined number of iterations N o = 10 6 , in each iteration, two standard normal variables are sampled independently with the intended number of samples N . The correlation between these two variables is then estimated and stored. This results in a vector of correlations of size N o that represents a distribution of empirically estimated correlations of independent variables and thus can be used to derive a cutoff based on some statistical measure. The 99.9 % percentile is considered in this work as the cutoff. As we are going to see later, N = 100 samples were used, and the calculated cutoff threshold was found to be 0.35 . This means that any correlation that, in absolute value, is larger than 0.35 is considered as statistically significant. The cross-sections that show at least one group to be statistically significant are analyzed at all their energy groups, as shown in the uncertainty quantification results section, in order to better understand the energy dependence of the estimated correlations.
  • The second aspect is the fact that if a very large number of inputs contribute approximately equally to the output, then it will be difficult to even distinguish them from the noise. This can be remedied by increasing the number of samples, since the cutoff threshold will decrease. Additionally, in most of the applications, few inputs contribute significantly more than the others, which facilitates their identification. The latter of course cannot be known in advance but only after performing the multi-physics calculations for the specific application.
  • The third and most important aspect is that the estimated correlations do not imply causation; in other terms, if an input shows strong correlation with an output, it does not mean that it is actually responsible for the output’s uncertainty. For this reason, a subjective analysis of the correlation matrices together with the physics understanding and experience can help deduce whether an input is actually important or it has high correlation, because it is correlated to another input, which is the one that is actually important. In any case, even if there is no objective way for such a simple uncertainty quantification approach to assign sensitivities in the multi-physics context, the resulting correlations are still very valuable and can inform future experiments to reduce the corresponding input uncertainties or correlations with other inputs. Finally, an ultimate test for this approach would be to compare the identified microscopic cross-sections for a steady-state standalone neutronics study with the ones identified using a perturbation-based approach. The validation of this proposed approach is left for future work.

3. LWR-UAM Phase III Core Study

3.1. Case Studies

The exercises I-2a, I-2b and I-2c of LWR-UAM Phase III [24] consist of various transient and steady-state calculations for the TMI-1 core. The TMI-1 core has a nominal power of 2772 MW th and consists of 177 UO 2 fuel assemblies with 11 different designs. The different designs consist of variations in UO 2 enrichment, number of gadolinium rods and number of discrete burnable poisons.
Exercise I-2a includes two steady-state calculations for BOC at HFP and Hot Zero Power (HZP) conditions. The main goal of this work is to apply and test the robustness of the multi-physics framework, so we focus on the HFP condition at both BOC and EOC. The 3D burnup distribution at BOC and EOC is provided in the benchmark specifications, and its axially averaged radial distribution is shown in Figure 10. The average core burnup is 18.08 GWD/MTU at BOC and 42.06 GWD/MTU at EOC. Exercise I-2b consists of a cycle depletion calculation from BOC to EOC. The cycle length is 664 Effective Full Power Days (EFPD). Exercise I-2c includes a series of REA from different initial conditions. The control rod with the highest worth, which is fully inserted into the core, is ejected at a constant rate of 2380 cm/s. No SCRAM is simulated to allow for a better understanding of the feedback effects. The location of the ejected control rod is the D12 in Figure 10. During the transient, the boron concentration and the position of the other control rods are assumed to be constant. Four different initial states are considered for the REA in the LWR-UAM benchmark: HZP BOC, HFP BOC, HZP EOC and HFP EOC. In this work, only the HFP BOC and EOC are modeled, because they lead to the most limiting condition in terms of the maximum fuel temperature and fuel stored enthalpy. The inserted reactivity in all cases is relatively small, and the HZP cases do not lead to large temperature variations.
The input uncertainty quantification for this work was performed based on the LWR-UAM Phase II [37] and Phase III [24] specifications and recommendations. It is important to mention that for the fuel and cladding thermal conductivities, the reduced uncertainty recommended in LWR-UAM Phase III and obtained after a Bayesian calibration is used in this work. A total of N = 100 macroscopic cross-sections samples (Section 3.2) for each unique assembly and reflector were generated with Sampler (Section 2.1.2). The Wilk’s formula was used to calculate that 93 samples were required to obtain a 95% tolerance interval with 95% confidence. A maximum of 1000 samples in SCALE can currently be generated from the SCALE libraries. Due to the computational cost for such uncertainty analysis, a sample size of N = 100 was deemed sufficient to test the computational framework. This sampling size choice was also justified in [38], where some of the LWR-UAM Phase III exercises were studied with a TRACE-PARCS coupling. More details about the Wilk’s formula and its limitations can be found in [39]. The sampled cross-sections were used together with five additional CTF thermal-hydraulics parameters, three boundary conditions and nine manufacturing parameters that are defined in Table 2 and Table 3 together with their probability density functions (PDF). In Table 2, the PDF are presented as multipliers on the nominal values, while in Table 3, the PDF are presented with the nominal values of the parameters along with their absolute standard deviations. The same standard deviation was used for all 235U enrichments that vary from 4.0 to 5.0 w/o%. The manufacturing parameters were originally sampled with Sampler using a random sampling size of N = 100 and were applied in the Polaris lattice calculations, while in CTF, a consistent perturbation of relevant impacted parameters was performed for each sample, as explained in Section 2.4. The perturbations for all manufacturing parameters were applied to all fuel rods in the assembly. The CTF parameters and the boundary conditions were sampled with Dakota, using a random sampling of size N = 100 .

3.2. Modeling

The core modeling is performed using a conventional two-step approach employed in LWR calculations. In the first step, lattice neutronics calculations are performed with Polaris. The ENDF/B-VII.1 cross-section library with SCALE’s 56-group structure is used to generate the two-group macroscopic cross-sections with a 0.625 eV energy cutoff. SCALE’s perturbation factors are used to generate the perturbed homoginized cross-section libraries. The perturbation factors are derived from the SCALE covariance library in 56 groups. The lattice calculations are performed for each unique fuel assembly design at 22 burnup points ranging from 0.0 to 65.0 GWD/MTU, 45 branch points and 8 history points. Fuel temperature, moderator density, boron concentration, and control insertion are the state parameters selected to define the branch and history points. The parameter coverage is given in Table 4 for the fuel assemblies. The four types of assembly layout are given in Figure 11.
Radial and axial reflector libraries were generated based on the LWR-UAM-Phase I specifications [40]. The radial reflector is constructed with three slabs neighboring a fuel assembly (Figure 12a). From left to right in Figure 12a is represented the fuel assembly, a 0.2 cm layer of coolant (first slab, blue), the shroud (second slab, yellow) and another layer of coolant (third slab, blue). The axial reflector is modeled with two slabs adjacent to a fuel assembly Figure 12b. From left to right in Figure 12b is represented the fuel assembly, the reflector (coolant, blue) and the core plate (second slab, yellow). The sub-regions in each slab represent the mesh grid used for the transport calculations in Polaris. The fuel assemblies in the radial and axial reflector models only provide a flux background using a critical spectrum to compute the homogenized cross-sections, so the reflector libraries only account for the reflector region. Reflective boundary conditions are implemented on the west side of the model, and vacuum boundary conditions are implemented on the east side of the model. The reflector libraries contain one history and 12 branches. The branch parameters for reflector libraries are the moderator density and boron concentration, and the same coverage as the fuel assemblies is used.
In the second step, the recently developed CTF-PARCS coupling was applied in an assembly-resolved model, leading to a total of 178 radial cells (177 fuel assemblies + 1 for the radial reflector). Axially, the core was divided into 28 meshes, 26 for the fuel active length and two for the top and bottom reflector. This level of nodalization is in line with LWR-UAM Phase III specifications and is visualized in Figure 13. In future uncertainty quantification studies, the uncertainties due to nodalization should also be considered using the Richardson extrapolation, as described in [41]. One representative fuel pin is modeled for each assembly with 10 radial nodes in the fuel and 2 radial nodes in the cladding. For each transient calculation, a steady-state coupled calculation is performed initially, and the core is brought to criticality by an eigenvalue fission source normalization. It is important to note that the normalization approach can have a non-negligible impact on the uncertainty quantification [42,43]. For the depletion calculation, a total of 24 depletion steps are performed. The first four steps are of 1 EFPD, and the next 20 steps are of 30 EFPD. The discussed modeling so far represents the standard fidelity of multi-physics calculations applied to uncertainty quantification studies, such as the one used in [22]. Compared to this standard fidelity and especially with regard to the thermal-hydraulics modeling in such multi-physics coupled calculations, some improvements in the thermal modeling have been implemented in this work, as discussed also in Section 2.3. The following improvements were implemented in CTF and used in the simulations presented in Section 4:
  • Consideration of 3D burnup distribution in the thermal calculations. The burnup impacts the thermal conductivity and thus the fuel temperature.
  • Radial power and burnup distribution within the representative fuel pin thermal calculation. The TUBRNP model, typically used in fuel performance codes, was used to calculate the burnup and power profiles for each local mesh based on the burnup of the mesh.
  • Exchange of the updated burnup distribution between PARCS and CTF in the depletion calculations.
  • Use of Santamarina effective Doppler temperature.

4. Results

In this section, we present the obtained results and statistical analysis for the different cases, using the previously discussed computational framework. The results are presented into three subsections. The first subsection is for the HFP BOC REA and includes two subdivisions: one for the steady-state results and one for the transient. The second subsection is for the HFP EOC REA and includes the same two subdivisions. The third and last subsection discusses the results for the cycle depletion from BOC to EOC. This creates a total of five distinctive case studies. For each case study, the results are analyzed in a hierarchy of functional complexity. At a first level, scalar quantities such as the maximum temperature during the transient are investigated with results about their Empirically estimated Cumulative Distributions (ECDF) and their first two statistical moments. For these quantities, the Spearman rank correlations are calculated, and as mentioned previously, only the correlations above the predefined absolute threshold of 0.35 are selected as statistically significant. It is important to stress that (a) these correlations do not necessarily imply causation, and (b) to derive any meaningful conclusion about sensitivities, a subjective judgment is necessary by understanding the underlying input uncertainties (e.g., correlations) and the physical phenomena modeled in each calculation. The second level of output quantities of interest are functional 1D quantities representing axial or temporal evolutions, such as the spatial maximum temperature trend during the transient, or the radially averaged axial power profile at HFP BOC. For these quantities, the mean is estimated together with the 95 % percentile as the upper bound and the 5 % percentile as the lower bound. The third level involves functional 2D radial output quantities such as the axially averaged power profile at HFP BOC. For these quantities, the mean and standard deviation are presented. All the selected scalar, functional 1D and functional 2D quantities for the steady state, transient and depletion cases are shown in Table 5, Table 6 and Table 7.

4.1. HFP BOC REA

4.1.1. Steady State

The steady-state calculation is performed for the HFP conditions as described in the specifications of LWR-UAM Phase III [24]. It is important to mention that an equilibrium Xenon and Samarium concentration is used. The results for the scalar outputs of interest are shown in Figure 14. The k eff shows an uncertainty of 0.44 % , which is a result in the same range of that found in the literature for typical LWR uncertainty quantification studies ( 0.5 % ) such as in [40]. Small uncertainties of 0.81 % and 0.37 % are obtained for P lin max and D cool min , respectively. The T fc max shows a larger uncertainty of 3.14 % .
In order to understand the sources of these uncertainties, we estimate the Spearman rank correlation coefficients. The 56-group cross-sections that exhibit strong correlations for the k eff and P lin max are shown in Figure 15. For the k eff , we see that the 235U ν ¯ p at low energies is the most significant input. This is in line with results from the LWR-UAM Phase I benchmark for steady state [40], since the ν ¯ p has a direct contribution to the produced neutrons and thus to the k eff . The positive correlation is reasonable, since an increase of the ν ¯ p will lead to an increase of the neutron population and then a higher k eff . As mentioned, correlation does not mean necessarily causation. This is the reason for example we cannot be sure that the ν ¯ p at all energies below 1 keV is important since they are highly correlated. We can, however, deduce that since among all the independent isotopic cross-sections only the ν ¯ p seems to be highly correlated and since that this is in accordance with our understanding of the physical phenomena, then we are quite confident that ν ¯ p is the most important input at low energies for the k eff . With a similar reasoning, we analyze the cross-section correlations for P lin max , where we can see a very strong correlation for the 238U elastic and inelastic scattering at high energies (above 1 MeV). From these two cross-sections, the inelastic scattering of 238U is a quantity that has been shown to be one of the dominant inputs in uncertainty quantification studies such as in the LWR-UAM Phase I benchmark [40]. For the elastic cross-section, we have little information, but as we can see in Figure 16, it is strongly negatively correlated with the inelastic scattering of 238U for this specific high-energy group. This is probably an indicator that the effect of both cross-sections is essentially one effect, but since they are highly correlated to each other, both appear to be statistically significant. To better understand the impact of these scattering cross-sections, in Figure 17, we show the estimated Spearman correlation between the axially maximum linear power in every assembly and the 238U elastic and inelastic cross-sections. In this figure, only the maximum correlation in absolute value is shown for each assembly. What can be observed is that the negative correlation is also evident in these radial results and that two distinctive regions are present. One region is in the center with positive/negative contributions for the elastic/inelastic cross-sections and one region is in the periphery with negative/positive contributions for the elastic/inelastic cross-sections. The highly negative correlations between the elastic and inelastic cross-sections do not allow a separation of their impact on the two regions. What can be deduced is that based on the location of the maximum power, the impact of these cross-sections can be reversed. This is indeed what we observe in all the next studies, because the maximum local power will move toward an assembly in the periphery of the core. We attribute the presence of these two regions to the fact that the flux is normalized to 1 and thus to the differential impact of the scattering on the distribution of neutrons in the core, leading essentially to more neutrons pushed toward the periphery. For both k eff and P lin max , no other source of uncertainty from the manufacturing, boundary conditions or CTF inputs was found to be statistically significant.
For the other outputs, no cross-section showed a statistically significant correlation, but instead for the T fc max , the following inputs and their respective Spearman correlations were found:
  • Fuel density: −0.39.
  • Cladding inner radius: 0.75.
  • Fuel thermal conductivity: −0.35.
  • Gap conductance: −0.45.
The most important input seems to be the cladding inner radius, which defines the gap thickness and thus gap conductance. A positive correlation is found because an increase of the gap width leads to a decrease of the gap conductance, less heat being extracted from the fuel resulting in a higher T fc max . For the same reason, the gap conductance is found to be the second significant input and with a negative correlation. Close to the cutoff threshold, we also find significant contributions from the fuel thermal conductivity and the fuel density, both with a negative sign. The former is simple to understand, because an increase in the thermal conductivity increases the heat extracted from the fuel and thus decreases the T fc max . For the latter, the effect is more difficult to explain, because it is an input to both Polaris and CTF. One explanation can be the contribution of the fuel density through the thermal conductivity, since an increase of the density will lead to an increase of the thermal conductivity. Another explanation could be the self-shielding on the fuel pellet outer layers and the effect on the neutron utilization: the change in fuel density in Polaris increases both the number density of 235U and 238U. The two isotopes have competitive effects on the neutron utilization, i.e., 235U will tend to increase the fission rate and 238U will increase the absorption rate due to resonances. The increase of 238U density may be dominant as compared to the 235U density. The two effects (conducitivity and neutronics) may be both contributing, or one of the two effects predominates. To verify those postulates, two tests can be performed for future work: (a) the density of the fuel can be modified in Polaris in a few samples, and the effect on T fc max can be evaluated without modifying the conductivity in the thermal-hydraulics calculations, and vice versa, (b) the fuel density in Polaris can be kept constant, but the fuel conductivity can be sampled in the thermal hydraulics calculations to evaluate the effect of the conductivity alone.
Regarding the D cool min , only the mass flux and the inlet temperature were found to be statistically significant. The impact of both is obvious for the D cool min and their Spearman correlations are the following:
  • Mass flux: 0.46.
  • Inlet temperature: −0.75.
The functional 1D results for the different axial profiles are shown in Figure 18. The mean is estimated together with the upper and lower bounds defined as the 95 % and 5 % percentiles, respectively. We can see that the linear power P lin , ax ave has a shape close to a cosine slightly shifted toward the bottom with a very small uncertainty. A similar trend is observed in [22]. The T fc , ax ave follows the same trend as P lin , ax ave with the upper and lower bounds ranging at ∼50 K . For the T cool , ax ave and D cool , ax ave , the expected trends are found with relatively small uncertainty bounds.
The functional 2D results are presented in Figure 19 and Figure 20 where the radial 2D mean and standard deviation of P lin , rad ave and T fc , rad ave are shown, respectively. We can see that the highest linear power is in the location E10 and all its symmetrical locations, while the largest uncertainty is ∼2.7% located in the central assembly. The centerline temperature profiles follow the linear power profile with the hot spot being at the location of E10. The highest uncertainty is ∼3.2% and is located in the fuel assemblies neighboring the central assembly.

4.1.2. Transient

The REA calculation is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24] starting from critical initial conditions reached by using fission source normalization. As for the steady state, an equilibrium Xenon and Samarium concentration is used. The results for the scalar outputs of interest are shown in Figure 21. The ρ i n has a mean value of 0.391$ and an uncertainty of 12.5 % . The P lin max mean has increased by a factor of 3 compared to the steady state with a significant larger uncertainty of 12.8 % . The T fc max follows the same behavior with an increase of ∼700 K over the steady-state average but with a smaller uncertainty of 1.7 % . This decrease is attributed to the fact that at steady state, the gap conductance could be open or closed depending on the initial gap width, while in the transient, the expansion of the fuel pin always leads to a gap closure and consequently decreases the impact of the gap width and the gap conductance. The H f max shows a mean value of 91 cal/g with an uncertainty of 2.0 % similar to the T fc max .
Investigating now the sources of uncertainties for these scalar output quantities, the Spearman rank correlation coefficients are calculated. Two clusters of trends can be observed coming from the cross-sections. The first cluster corresponds to ρ i n and P lin max where the 56-group cross-sections show the same strong correlations for both of them as can be seen in Figure 22a,b. Similar to the steady-state results, the 238U elastic and inelastic scattering cross-sections at high energies are statistically significant. Interestingly enough, they show an inverse sign for the correlation with regard to the steady-state results. This is attributed to the different location of the maximum linear power between the steady-state and transient results. The location moves from the center toward the periphery, and as can be seen in Figure 17, this leads to an inverse correlation. Although these cross-sections show strong correlations, the 238U ν ¯ d and ν ¯ p at high energies show even stronger correlations. The opposing signs between the ν ¯ d and ν ¯ p are attributed to their strong negative correlation in these energies, as can be seen in Figure 23. Further analysis of the 238U ν ¯ d and ν ¯ p uncertainties leads to the conclusion that the effective delayed neutron fraction is strongly correlated with these cross-sections and specifically with opposing correlations signs compared to ρ i n . This means that probably the 238U ν ¯ d and ν ¯ p cross-sections determine the uncertainty of the effective delayed neutron fraction, which is in the denominator of the ρ i n . What we can deduce at the end from this analysis is that the uncertainty of the ρ i n and P lin max is determined by the ejected control rod worth and by the effective delayed neutron fraction. The former uncertainty is mainly determined by 238U inelastic and elastic scattering uncertainties, while the latter is mainly determined by the 238U ν ¯ d and ν ¯ p uncertainties. The remaining two quantities T fc max and H f max belong to the second cluster of trends and the cross-sections with stronger correlations are only the 238U elastic and inelastic scattering cross-sections, as can be seen in Figure 22c,d.
From all the other sources of uncertainties, only the fuel density is statistically significant for the T fc max and H f max with a Spearman rank correlation of −0.37 for both. As explained in the steady-state calculations, this correlation is attributed to the impact of the fuel density on the thermal conductivity and/or the absorption rate of 238U. The fact that the gap conductance and gap width do not seem to be significant reinforces the explanation of the reduced uncertainties on the T fc max compared to the steady state.
The functional 1D results over time for P lin , t max , T fc , t max and H f , t max are presented in Figure 24. It can be seen that P lin , t max shows large uncertainties around the maximum, as seen also in the scalar results, but then the uncertainty decreases rapidly. The location of the maximum linear power is located during the whole transient in the assemblies D13 and E12 that are neighboring the control rod ejection assembly (D12). The T fc , t max behavior shows that the maximum fuel centerline temperature stabilizes after 35 s . Initially, the hot spot location coincides with the location of the maximum linear power, but at around 10 s , as indicated by the non-linearity in the behavior, the hot spot location moves to the assembly where the rod was ejected, D12. This happens because the D12 assembly has higher burnup, as can be seen in Figure 10, and thus, the deteriorating thermal conductivity leads to larger temperatures later on in the transient. The uncertainty of T fc , t max creates bounds of ∼80 K . The H f , t max behavior follows closely the behavior of T fc , t max .
The functional 2D results at the end of the transient for P lin , end ave and T fc , end ave are presented in Figure 25 and Figure 26. The P lin , end ave shows the expected peaking around the location of the control ejection. The same trend is observed for T fc , end ave . Concerning the uncertainties, the P lin , end ave shows a standard deviation up to 3.0 % in the central assembly, while the T fc , end ave shows a standard deviation of up to 3.5 % in the locations of the fresh assemblies with 0 GWD/MTU burnup.

4.2. HFP EOC REA

4.2.1. Steady State

The steady-state calculation is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24] at the EOC. The results for the scalar outputs of interest are shown in Figure 27. The k eff shows a slight increase compared to the BOC case, while the P lin max shows a significant increase. The uncertainty of D cool min remains similar, while the T fc max decreases almost by half.
The underlying sources of these uncertainties are investigated through the estimation of the Spearman rank correlation coefficients. Some of the 56-group cross-sections show strong correlations for the k eff , P lin max and T fc max as shown in Figure 28, while no cross-section is identified as significant for D cool min . For k eff , a very interesting observation can be made since we do not see anymore the 235U ν ¯ p being the dominant input, but instead, we see contributions from the 238U elastic and inelastic scattering at high energies and 239Pu fission and capture at low energies. This is attributed to the increased average burnup of the core and thus to stronger 239Pu buildup that contributes a significant part of the neutrons production. This explains also the positive sign for the fission cross-section and the negative sign for the capture cross-section. The results for P lin max are similar to the BOC case with the 238U elastic and inelastic scattering at high energies being the statistically significant inputs. The only difference is the inverse sign of the correlations, which is something attributed to the same reason as for the BOC transient results. The location of the maximum moves from the center region (E10) at BOC to the periphery (C10) at EOC. These regions, as shown in Figure 17, have opposite correlations for the 238U elastic and inelastic scattering. For T fc max , we see the same cross-sections as important, but in the BOC calculations, these were not identified. This is because probably at BOC, their correlation was below the cutoff threshold, while at EOC, it increased slightly above the threshold.
From the other inputs, none was found statistically significant for the k eff and P lin max . For T fc max , the following were found to be significant:
  • Fuel density: −0.45.
  • Fuel thermal conductivity: −0.46
The fuel density and thermal conductivity correlation increased for T fc max compared to BOC, while we do not see as important anymore the gap heat transfer-related inputs. This is attributed to the increased average burnup at the EOC. Regarding the D cool min , the same inputs with the BOC case are found significant with their correlation being:
  • Mass flow rate: 0.41.
  • Inlet temperature: −0.79.
The functional 1D results for the different axial profiles are shown in Figure 29. The main impact is on the behavior of P lin , ax ave and T fc , ax ave , where a strong bottom peaking is observed. A similar trend is observed in [22]. The uncertainties on all the functional quantities are small, as can be seen from the tight upper and lower bounds in Figure 29.
The functional 2D results are presented in Figure 30 and Figure 31 where the radial 2D mean and standard deviation of P lin , rad ave and T fc , rad ave are shown, respectively. We can see that the peaking of the linear power is in the location E10 similar to the BOC case but also in some peripheral assemblies such as C10 assembly. The P lin , rad ave uncertainty is up to ∼2.5% and in the central assembly. The T fc , rad ave follows the mean profile of P lin , rad ave and exhibits uncertainties up to ∼1.4% located in the fuel assemblies in the periphery of the core.

4.2.2. Transient

The REA calculation is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24] at EOC and starting from critical initial conditions using fission source normalization. The results for the scalar outputs of interest are shown in Figure 32. The ρ i n mean value has decreased to 0.357% but its uncertainty increased to 16.4 % . The P lin max mean has decreased compared to the BOC REA due to the smaller inserted reactivity, and its uncertainty remained similar at 11.9 % . The T fc max shows a significant decrease of ∼400 K compared to the BOC REA with a similar uncertainty of 1.7 % . This decrease is attributed to a less extreme REA transient at EOC conditions. The H f max decreases significantly its mean to 76 cal / g following the T fc max trend with a similar uncertainty of 1.6 % .
In order to understand the change in the obtained uncertainties, the Spearman rank correlation coefficients are calculated. Similar to the BOC REA, two clusters of trends can be observed concerning the cross-sections. The first cluster corresponds to ρ i n and P lin max , where the 56-group cross-sections that exhibit strong correlations are shown in Figure 22a,b. The difference with respect to BOC REA is that the 238U elastic and inelastic cross-sections at high energies contributions have decreased significantly, while the 239Pu ν ¯ d at low energies shows a strong statistical significance. Apart from that, the dominant cross-sections are still the 238U ν ¯ d and ν ¯ p at high energies. For T fc max and H f max , the same cross-sections as for the BOC REA case appear, namely, the 238U elastic and inelastic cross-section at high energies, as can be seen in Figure 33c,d.
From all the other sources of uncertainties, only the fuel density is statistically significant for the T fc max with a Spearman rank correlation of −0.35. As explained in the HFP BOC case, this correlation is attributed to the impact of the fuel density on the thermal conductivity.
The functional 1D results over time for P lin , t max , T fc , t max and H f , t max are presented in Figure 34, and similar conclusions to the BOC REA can be reached. The P lin , t max shows large uncertainties around the maximum with the uncertainty decreasing rapidly afterwards. The location of the maximum linear power is located during the whole transient in the assemblies D13 and E12 that are neighboring the control rod ejection assembly (D12). The T fc , t max behavior shows that the maximum fuel centerline temperature stabilizes after 30 s and is now located at the assembly E12. The uncertainty of T fc , t max creates bounds of ∼40 K . The H f , t max behavior follows closely the behavior of T fc , t max .
The functional 2D results at the end of the transient for P lin , end ave and T fc , end ave are presented in Figure 35 and Figure 36. The P lin , end ave peaking moves toward the center of the core compared to the BOC REA. This applies to the T fc , end ave as well that now follows more closely the P lin , end ave behavior. Concerning the uncertainties, the P lin , end ave shows a similar standard deviation up to 2.5 % in the central assembly, while the T fc , end ave decreased by half a standard deviation of up to 1.5 % in both central and peripheral assemblies.

4.3. Depletion

The depletion from BOC to EOC calculations is performed for the HFP conditions described in the specifications of LWR-UAM Phase III [24]. The results for the scalar outputs of interest are shown in Figure 37. The obtained Bu ave mean is 40.381 GWD/MTU, which is a value relatively close to the reference value of 42.06 GWD/MTU found in the LWR-UAM specifications. The Bu max mean is 62.78 GWD/MTU. The uncertainties for Bu ave and Bu max are 0.3 % and 0.5 % , respectively.
The Spearman rank correlation coefficients between the 56-group cross-sections and Bu ave do not show any statistical significant cross-section. For the Bu max , however, four cross-sections are found significant and are shown in Figure 38. The main cross-sections are the 238U elastic and inelastic scattering at high energies. The sign of these correlations is similar to the BOC steady-state results and reversed compared to the BOC transient and EOC steady state and transient. This, as explained in Figure 17, is due to the location of the maximum that occurs in the central region of the core. At a second order, the 235U elastic and fission at low energies are identified, but they are very close to the cutoff threshold. From the remaining inputs, only the fuel density is significant, with a dominant role for Bu ave showing a correlation of −0.91, while for the Bu max , it has a correlation of −0.41. The impact of the fuel density and its negative sign is obvious, because the burnup step calculated by PARCS has the fuel density in the denominator.
The functional 1D results are split into axial and temporal. The results for the temporal quantities are shown in Figure 39. As we can see, the Bu t max is increasing linearly for most of the depletion with small uncertainty bounds. The P lin , t max decreases during the depletion with small uncertainties indicating the decrease of the axial offset. The T fc , t max evolution over depletion indicates uncertain bounds of ∼80 K . The T fc , t max trend shows a decrease from 1700 K to 1540 K during half of the depletion and then a stabilization around this value. The results for the axial quantities are shown in Figure 40. We can see that the linear power P lin , BOC ave has a shape close to a cosinus with a slight shift toward the bottom and exhibits a very small uncertainty. At the EOC, the P lin , EOC ave shows a double peaking profile with one peak on the bottom and one on the top with very small uncertainties. The obtained radially averaged axial burnup profile Bu ax , EOC ave is compared in Figure 40c with the reference results of the LWR-UAM specifications. The comparison indicates a relatively good agreement for most of the fuel active height and larger discrepancies close to the top and bottom reflectors.
The functional 2D radial results for the obtained axially averaged burnup distribution Bu rad , EOC ave are presented in Figure 41. The estimated standard deviation reaches values up to 1.2 % in locations close to the periphery. In Figure 42, the discrepancy between the nominal Bu rad , EOC ave and the reference from the LWR-UAM specifications is shown. The errors span from 3.5 % to 1.5 % , indicating a reasonable agreement.

5. Conclusions

In this work, the developed CTF-PARCS core multi-physics computational framework for efficient Light Water Reactor (LWR) steady-state, depletion and transient uncertainty quantification is presented. Novel features of this framework include a versatile Python pre-processing script for streamlining the macroscopic cross-section generation process and the capability to estimate correlations between outputs of interest and the 56-group microscopic cross-sections, fission yields and all other multi-physics input parameters.
The developed computational framework is applied to three exercises of LWR-UAM Phase III that cover steady-state, depletion and transient calculations. The results show that the maximum fuel centerline temperature upper bound, defined as the 95 % percentile, during the Beginning Of Cycle (BOC) Rod Ejection Accident (REA) is ∼2600 K , while during the End Of Cycle (EOC) REA, it is ∼2200 K . The results, in terms of correlations with the 56-group microscopic cross-sections, show that few isotope-reactions are found as statistically significant with strong correlations:
  • At BOC steady state, the 235U ν ¯ p at low energies (<1 keV) and the 238U inelastic and elastic scattering at high energies (>1 MeV, which corresponds to the 238U inelastic scattering threshold).
  • At BOC REA, the 238U inelastic and elastic scattering at high energies as well as the 238U ν ¯ p and ν ¯ d at high energies.
  • At EOC steady state, the 238U inelastic and elastic scattering at high energies as well as 239Pu fission and capture cross-section at low energies (<10 eV).
  • At EOC REA, the 238U inelastic and elastic scattering at high energies, the 238U ν ¯ p and ν ¯ d , and 239Pu ν ¯ d at low energies (<100 eV).
The results highlight the efficiency and robustness of the computational framework and the ability to extract information, although subjective, up to the microscopic cross-sections, which is something that can be very useful for identifying future experiments to reduce the nuclear data uncertainties.
Future steps based on this work will involve validating the computational framework and estimating the impact of different modeling approaches. Furthermore, parallel domain decomposition could be implemented in CTF to accelerate the calculations and allow even a pin-cell level of multi-physics coupling. The uncertainty quantification approach could also be investigated more thoroughly in order to derive a more objective way of discussing the obtained correlations for the microscopic cross-sections.

Author Contributions

Conceptualization, G.K.D., P.R. and A.A.; methodology, G.K.D., P.R., A.A. and K.I.; software, G.K.D., P.R., A.A. and J.H.; validation, P.R., A.A., K.I. and M.A.; formal analysis, K.I., M.A. and J.H.; investigation, G.K.D. and P.R.; resources, K.I., M.A. and J.H.; data curation, G.K.D., P.R. and A.A.; writing—original draft preparation, G.K.D., P.R. and A.A.; writing—review and editing, K.I., M.A. and J.H.; visualization, G.K.D., P.R. and A.A.; supervision, K.I., M.A. and J.H.; project administration, K.I., M.A. and J.H.; funding acquisition, K.I., M.A. and J.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
BEPUBest Estimate Plus Uncertainty
BWRBoiling Water Reactor
BOCBeginning of Cycle
CASL   Consortium for Advanced Simulation of Light water reactors
DOEDepartment of Energy
E.C.European Commission
EFPDEffective Full Power Days
EOCEnd of Cycle
HFPHot Full Power
HPCHigh Performance Computing
HZPHot Zero Power
LWRLight Water Reactor
MPIMessage Passing Interface
NEANuclear Energy Angency
NCSUNorth Carolina State University
OECDOrganisation for Economic Co-operation and Development
PWRPressurized Water Reactor
REARod Ejection Accident
SCRAMSafety Control Rod Axe Man
TMIThree Mile Island
U.S.United States
VVERWater Water Energetic Reactor

References

  1. Rohatgi, U.S.; Kaizer, J.S. Historical perspectives of BEPU research in US. Nucl. Eng. Des. 2020, 358, 110430. [Google Scholar] [CrossRef]
  2. Kochunas, B.; Collins, B.; Stimpson, S.; Salko, R.; Jabaay, D.; Graham, A.; Liu, Y.; Seog Kim, K.; Wieselquist, W.; Godfrey, A.; et al. VERA Core Simulator Methodology for Pressurized Water Reactor Cycle Depletion. Nucl. Sci. Eng. 2017, 185, 217–231. [Google Scholar] [CrossRef]
  3. Lefebvre, R.A.; Langley, B.R.; Miller, P.; Delchini, M.; Baird, M.L.; Lefebvre, J.P. NEAMS Workbench Status and Capabilities; Technical Report ORNL/TM-2019/1314; Oak Ridge National Laboratory, Reactor Nuclear Systems Division: Oak Ridge, TN, USA, 2019. [Google Scholar]
  4. Permann, C.J.; Gaston, D.R.; Andrš, D.; Carlsen, R.W.; Kong, F.; Lindsay, A.D.; Miller, J.M.; Peterson, J.W.; Slaughter, A.E.; Stogner, R.H.; et al. MOOSE: Enabling massively parallel multiphysics simulation. SoftwareX 2020, 11, 100430. [Google Scholar] [CrossRef]
  5. Targa, A. Development of Multi-Physics and Multi-Scale Best Effort Modelling of Pressurized Water Reactor under Accidental Situations. Ph.D. Thesis, Université Paris Saclay (COmUE), Saclay, France, 2017. [Google Scholar]
  6. Fiorina, C.; Clifford, I.; Aufiero, M.; Mikityuk, K. GeN-Foam: A novel OpenFOAM® based multi-physics solver for 2D/3D transient analysis of nuclear reactors. Nucl. Eng. Des. 2015, 294, 24–37. [Google Scholar] [CrossRef]
  7. Valtavirta, V.; Hovi, V.; Loukusa, H.; Rintala, A.; Sahlberg, V.; Tuominen, R.; Leppänen, J. Kraken: An upcoming Finnish reactor analysis framework. In Proceedings of the International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2019), 25–29 August 2019; American Nuclear Society (ANS): Portland, OR, USA, 2019; pp. 786–795. [Google Scholar]
  8. Cherezov, A.; Park, J.; Kim, H.; Choe, J.; Lee, D. A Multi-Physics Adaptive Time Step Coupling Algorithm for Light-Water Reactor Core Transient and Accident Simulation. Energies 2020, 13, 6374. [Google Scholar] [CrossRef]
  9. Hou, J.; Avramova, M.; Ivanov, K. Best-Estimate Plus Uncertainty Framework for Multiscale, Multiphysics Light Water Reactor Core Analysis. Sci. Technol. Nucl. Install. 2020, 2020, 7526864. [Google Scholar] [CrossRef]
  10. Salko, R., Jr.; Avramova, M.; Wysocki, A.; Toptan, A.; Hu, J.; Porter, N.; Blyth, T.S.; Dances, C.A.; Gomez, A.; Jernigan, C.; et al. CTF 4.0 Theory Manual; CASL, Oak Ridge National Lab: Oak Ridge, TN, USA, 2019. [Google Scholar] [CrossRef]
  11. Downar, T.; Barber, D.; Matthew Miller, R.; Lee, C.; Kozlowski, T.; Lee, D.; Xu, Y.; Gan, J.; Joo, H.; Cho, J.; et al. PARCS: Purdue advanced reactor core simulator. In Proceedings of the PHYSOR 2002—International Conference on the New Frontiers of Nuclear Technology, 7–10 October, Seoul, Republic of Korea; American Nuclear Society: La Grange Park, IL, USA, 2002. [Google Scholar]
  12. Wang, J.; Wang, Q.; Ding, M. Review on neutronic/thermal-hydraulic coupling simulation methods for nuclear reactor analysis. Ann. Nucl. Energy 2020, 137, 107165. [Google Scholar] [CrossRef]
  13. Ivanov, K.; Avramova, M. Challenges in coupled thermal–hydraulics and neutronics simulations for LWR safety analysis. Ann. Nucl. Energy 2007, 34, 501–513. [Google Scholar] [CrossRef]
  14. Ragusa, J.C.; Mahadevan, V.S. Consistent and accurate schemes for coupled neutronics thermal-hydraulics reactor analysis. Nucl. Eng. Des. 2009, 239, 566–579. [Google Scholar] [CrossRef]
  15. Senecal, J.P.; Ji, W. Development of an efficient tightly coupled method for multiphysics reactor transient analysis. Prog. Nucl. Energy 2018, 103, 33–44. [Google Scholar] [CrossRef]
  16. Zhang, H.; Guo, J.; Lu, J.; Niu, J.; Li, F. A comparison of coupling algorithms for N/TH transient problems in HTR. In Proceedings of the International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2017), Jeju, Korea, 16–20 April 2017. [Google Scholar]
  17. Zhang, H.; Guo, J.; Lu, J.; Li, F.; Xu, Y.; Downar, T.J. An Assessment of Coupling Algorithms in HTR Simulator TINTE. Nucl. Sci. Eng. 2018, 190, 287–309. [Google Scholar] [CrossRef]
  18. Hamilton, S.; Berrill, M.; Clarno, K.; Pawlowski, R.; Toth, A.; Kelley, C.; Evans, T.; Philip, B. An assessment of coupling algorithms for nuclear reactor core physics simulations. J. Comput. Phys. 2016, 311, 241–257. [Google Scholar] [CrossRef] [Green Version]
  19. Ellis, M.; Watson, J.; Ivanov, K. Progress in the development of an implicit steady state solution in the coupled code TRACE/PARCS. Prog. Nucl. Energy 2013, 66, 1–12. [Google Scholar] [CrossRef]
  20. Ivanov, K.; Avramova, M.; Kamerow, S.; Kodeli, I.; Sartori, E.; Ivanov, E.; Cabellos, O. Benchmarks for Uncertainty Analysis in Modelling (UAM) for the Design, Operation and Safety Analysis of LWRs, Volume I: Specification and Support Data for Neutronics Cases (Phase I); OECD Nuclear Energy Agency: Paris, France, 2013. [Google Scholar]
  21. Delipei, G.K. Development of an Uncertainty Quantification methodology for Multi-Physics Best Estimate Analysis and Application to the Rod Ejection Accident in a Pressurized Water Reactor. Ph.D. Thesis, Université Paris Saclay (COmUE), Saclay, France, 2019. [Google Scholar]
  22. Kaiyue, Z. Uncertainty Analysis Framework for the Multi-Physics Light Water Reactor Simulation. Ph.D. Thesis, North Carolina State University, Raleigh, NC, USA, 2020. [Google Scholar]
  23. U.S. NRC. TRACE V5.0 Theory Manual, Field Equations, Solution Methods, and Physical Models; Technical Report; U.S. NRC: Washington, DC, USA, 2007. [Google Scholar]
  24. Hou, J.; Avramova, M.; Ivanov, K.; Royer, E.; Jessee, M.; Zhang, J.; Wieselquist, W.; Pasichnyk, I.; Zwermann, W.; Velkov, K.; et al. Benchmarks for Uncertainty Analysis in Modelling (UAM) for the Design, Operation and Safety Analysis of LWRs, Volume III: Specification and Support Data for the System Cases (Phase III); OECD Nuclear Energy Agency: Paris, France, 2019. [Google Scholar]
  25. Jessee, M.A.; Wieselquist, W.A.; Evans, T.M.; Hamilton, S.P.; Jarrell, J.J.; Kim, K.S.; Lefebvre, J.P.; Lefebvre, R.A.; Mertyurek, U.; Thompson, A.B.; et al. POLARIS: A New Two-Dimensional Lattice Physics Analysis Capability for the SCALE Code System; Oak Ridge National Laboratory (ORNL): Oak Ridge, TN, USA, 2014. [Google Scholar]
  26. Rearden, B.T.; Jessee, M.A. SCALE Code System, Version 6.2.3; Technical Report ORNL/TM-2005/39; Oak Ridge National Laboratory, Reactor and Nuclear Systems Division: Oak Ridge, TN, USA, 2018. [Google Scholar]
  27. Chadwick, M.B.; Obložinský, P.; Dunn, M.; Danon, Y.; Kahler, A.; Smith, D.; Pritychenko, B.; Arbanas, G.; Arcilla, R.; Brewer, R.; et al. ENDF/B-VII.1 Nuclear Data for Science and Technology: Cross Sections, Covariances, Fission Product Yields and Decay Data. Nucl. Data Sheets 2011, 112, 2887–2996. [Google Scholar] [CrossRef]
  28. Ward, A.; Xu, Y.; Downar, T. GenPMAXS-v6.1.3—Code for Generating the PARCS Cross Section Interface File PMAXS; Technical Report; University of Michigan: Ann Arbor, MI, USA, 2015. [Google Scholar]
  29. U.S. DOE. CASL Phase II Summary Report; Technical Report CASL-U-2020-1974-000; ORNL/SPR-2020/1759; CASL; 2020. Available online: https://casl.gov/wp-content/uploads/2020/11/CASL_FINAL_REPORT_09.30.2020-002.pdf (accessed on 13 July 2022).
  30. European Commission. NURESAFE Project Final Report. Technical Report GA N° 323263; E.C.. 2016. Available online: https://cordis.europa.eu/docs/results/323/323263/final1-nuresafe-final-report-publishable-summary-june-29-2016.pdf (accessed on 13 July 2022).
  31. Dalbey, K.; Eldred, M.S.; Geraci, G.; Jakeman, J.D.; Maupin, K.A.; Monschke, J.A.; Seidl, D.T.; Swiler, L.P.; Tran, A.; Menhorn, F.; et al. Dakota A Multilevel Parallel Object-Oriented Framework for Design Optimization Parameter Estimation Uncertainty Quantification and Sensitivity Analysis: Version 6.12 Theory Manual; Sandia National Laboratories: Albuquerque, NM, USA, 2020. [Google Scholar] [CrossRef]
  32. Abarca, A.; Avramova, M.; Kostadin, I. Development of a General MPI Coupling Interface for Multi-Physics Analysis; ANS International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering: Lucca, Italy, 2021. [Google Scholar]
  33. Rowlands, G. Resonance absorption and non-uniform temperature distributions. J. Nucl. Energy. Parts A/B. React. Sci. Technol. 1962, 16, 235–236. [Google Scholar] [CrossRef]
  34. Bernard, D.; Calame, A.; Palau, J.M. LWR-UOx doppler reactivity coefficient best estimate plus (nuclear and atomic sources of) uncertainties. In Proceedings of the M&C—2017 International Conference on Mathematics & Computational Methods Applied to Nuclear Science & Engineering, Jeju, Korea, 16–20 April 2017. [Google Scholar]
  35. Lassmann, K.; O’Carroll, C.; van de Laar, J.; Walker, C. The radial distribution of plutonium in high burnup UO2 fuels. J. Nucl. Mater. 1994, 208, 223–231. [Google Scholar] [CrossRef]
  36. Iooss, B.; Marrel, A. An efficient methodology for the analysis and modeling of computer experiments with large number of inputs. In Proceedings of the UNCECOMP 2017 2nd ECCOMAS Thematic Conference on Uncertainty Quantification in Computational Sciences and Engineering, Rhodes Island, Greece, 15–17 June 2017; pp. 187–197. [Google Scholar]
  37. Hou, J.; Blyth, T.; Porter, N.; Avramova, M.; Ivanov, K.; Royer, E.; Sartori, E.; Cabellos, O.; Feroukhi, H.; Ivanov, E. Benchmarks for Uncertainty Analysis in Modelling (UAM) for the Design, Operation and Safety Analysis of LWRs, Volume II: Specification and Support Data for the Core Cases (Phase II); OECD Nuclear Energy Agency: Paris, France, 2019. [Google Scholar]
  38. Zeng, K.; Hou, J.; Ivanov, K.; Jessee, M.A. Uncertainty Quantification and Propagation of Multiphysics Simulation of the Pressurized Water Reactor Core. Nucl. Technol. 2019, 205, 1618–1637. [Google Scholar] [CrossRef]
  39. Porter, N.W. Wilks’ formula applied to computational tools: A practical discussion and verification. Ann. Nucl. Energy 2019, 133, 129–137. [Google Scholar] [CrossRef]
  40. Delipei, G.K.; Hou, J.; Avramova, M.; Rouxelin, P.; Ivanov, K. Summary of comparative analysis and conclusions from OECD/NEA LWR-UAM benchmark Phase I. Nucl. Eng. Des. 2021, 384, 111474. [Google Scholar] [CrossRef]
  41. Roy, C.J.; Oberkampf, W.L. A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing. Comput. Methods Appl. Mech. Eng. 2011, 200, 2131–2144. [Google Scholar] [CrossRef]
  42. Dokhane, A.; Vasiliev, A.; Hursin, M.; Rochman, D.; Ferroukhi, H. A critical study on best methodology to perform UQ for RIA transients and application to SPERT-III experiments. Nucl. Eng. Technol. 2021, 54, 1804–1812. [Google Scholar] [CrossRef]
  43. Delipei, G.K.; Garnier, J.; Le Pallec, J.; Normand, B. Multi-Physics Uncertainties Propagation in a PWR Rod Ejection Accident Modeling—Analysis Methodology and First Results; ANS Best Estimate Plus Uncertainty International Conference (BEPU 2018): Lucca, Italy, 2018. [Google Scholar]
Figure 1. Construction of the Polaris inputs from template.
Figure 1. Construction of the Polaris inputs from template.
Energies 15 05226 g001
Figure 2. On-the-fly verification of Polaris library generation.
Figure 2. On-the-fly verification of Polaris library generation.
Energies 15 05226 g002
Figure 3. General scheme of the MPI driven multi-physics multi-scale general simulation framework developed at NCSU.
Figure 3. General scheme of the MPI driven multi-physics multi-scale general simulation framework developed at NCSU.
Energies 15 05226 g003
Figure 4. CTF-PARCS external coupling using a client/server MPI Coupling Interface.
Figure 4. CTF-PARCS external coupling using a client/server MPI Coupling Interface.
Energies 15 05226 g004
Figure 5. Coupled code information for (a) feedback parameters and (a) transient temporal scheme.
Figure 5. Coupled code information for (a) feedback parameters and (a) transient temporal scheme.
Energies 15 05226 g005
Figure 6. Details of the coupled steady-state convergence.
Figure 6. Details of the coupled steady-state convergence.
Energies 15 05226 g006
Figure 7. Scheme of the coupled code convergence during depletion simulations.
Figure 7. Scheme of the coupled code convergence during depletion simulations.
Energies 15 05226 g007
Figure 8. Details of the coupled depletion convergence.
Figure 8. Details of the coupled depletion convergence.
Energies 15 05226 g008
Figure 9. Uncertainty quantification approach.
Figure 9. Uncertainty quantification approach.
Energies 15 05226 g009
Figure 10. Axially averaged 2D burnup distribution at (a) BOC and (b) EOC.
Figure 10. Axially averaged 2D burnup distribution at (a) BOC and (b) EOC.
Energies 15 05226 g010
Figure 11. 1/4 Assembly layouts modeled in Polaris. (a) Assembly without burnable poisons or UO 2 - Gd 2 O 3 rods; (b) Assembly with 4 UO 2 - Gd 2 O 3 rods and without burnable poisons; (c) Assembly with 8 UO 2 - Gd 2 O 3 rods and without burnable poisons; (d) Assembly with 8 UO 2 - Gd 2 O 3 rods and with burnable poisons.
Figure 11. 1/4 Assembly layouts modeled in Polaris. (a) Assembly without burnable poisons or UO 2 - Gd 2 O 3 rods; (b) Assembly with 4 UO 2 - Gd 2 O 3 rods and without burnable poisons; (c) Assembly with 8 UO 2 - Gd 2 O 3 rods and without burnable poisons; (d) Assembly with 8 UO 2 - Gd 2 O 3 rods and with burnable poisons.
Energies 15 05226 g011
Figure 12. Reflector configurations in Polaris, (a) radial reflector, (b) axial reflector.
Figure 12. Reflector configurations in Polaris, (a) radial reflector, (b) axial reflector.
Energies 15 05226 g012
Figure 13. Radial (a) and axial (b) nodalization of the CTF-PARCS modeling.
Figure 13. Radial (a) and axial (b) nodalization of the CTF-PARCS modeling.
Energies 15 05226 g013
Figure 14. ECDF of (a) k eff , (b) P lin max , (c) T fc max and (d) D cool min .
Figure 14. ECDF of (a) k eff , (b) P lin max , (c) T fc max and (d) D cool min .
Energies 15 05226 g014
Figure 15. Statistically significant cross-sections Spearman correlation coefficients for (a) k eff and (b) P lin max .
Figure 15. Statistically significant cross-sections Spearman correlation coefficients for (a) k eff and (b) P lin max .
Energies 15 05226 g015
Figure 16. SCALE 56-group correlation matrix between 238U elastic (MT = 2) and inelastic (MT = 4) scattering.
Figure 16. SCALE 56-group correlation matrix between 238U elastic (MT = 2) and inelastic (MT = 4) scattering.
Energies 15 05226 g016
Figure 17. Radial distribution of the maximum Spearman correlation coefficient between the axially maximum linear power and (a) 238U elastic scattering and (b) 238U inelastic scattering.
Figure 17. Radial distribution of the maximum Spearman correlation coefficient between the axially maximum linear power and (a) 238U elastic scattering and (b) 238U inelastic scattering.
Energies 15 05226 g017
Figure 18. Axial profile for (a) P lin , ax ave , (b) T fc , ax ave , (c) T cool , ax ave and (d) D cool , ax ave .
Figure 18. Axial profile for (a) P lin , ax ave , (b) T fc , ax ave , (c) T cool , ax ave and (d) D cool , ax ave .
Energies 15 05226 g018
Figure 19. Radial profile of P lin , rad ave (a) mean and (b) standard deviation.
Figure 19. Radial profile of P lin , rad ave (a) mean and (b) standard deviation.
Energies 15 05226 g019
Figure 20. Radial profile of T fc , rad ave (a) mean and (b) standard deviation.
Figure 20. Radial profile of T fc , rad ave (a) mean and (b) standard deviation.
Energies 15 05226 g020
Figure 21. ECDF of (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Figure 21. ECDF of (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Energies 15 05226 g021
Figure 22. Statistically significant cross-sections Spearman correlation coefficients for (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Figure 22. Statistically significant cross-sections Spearman correlation coefficients for (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Energies 15 05226 g022
Figure 23. SCALE 56-group correlation matrix between 238U ν ¯ p (MT = 456) and ν ¯ d (MT = 455).
Figure 23. SCALE 56-group correlation matrix between 238U ν ¯ p (MT = 456) and ν ¯ d (MT = 455).
Energies 15 05226 g023
Figure 24. Temporal evolution of the (a) P lin , t max , (b) T fc , t max and (c) H f , t max .
Figure 24. Temporal evolution of the (a) P lin , t max , (b) T fc , t max and (c) H f , t max .
Energies 15 05226 g024
Figure 25. Radial profile of P lin , end ave (a) mean and (b) standard deviation.
Figure 25. Radial profile of P lin , end ave (a) mean and (b) standard deviation.
Energies 15 05226 g025
Figure 26. Radial profile of T fc , end ave (a) mean and (b) standard deviation.
Figure 26. Radial profile of T fc , end ave (a) mean and (b) standard deviation.
Energies 15 05226 g026
Figure 27. ECDF of (a) k eff , (b) P lin max , (c) T fc max and (d) D cool min .
Figure 27. ECDF of (a) k eff , (b) P lin max , (c) T fc max and (d) D cool min .
Energies 15 05226 g027
Figure 28. Statistically significant cross-sections Spearman correlation coefficients for (a) k eff , (b) P lin max and (c) T fc max .
Figure 28. Statistically significant cross-sections Spearman correlation coefficients for (a) k eff , (b) P lin max and (c) T fc max .
Energies 15 05226 g028
Figure 29. Axial profile for (a) P lin , ax ave , (b) T fc , ax ave , (c) T cool , ax ave and (d) D cool , ax ave .
Figure 29. Axial profile for (a) P lin , ax ave , (b) T fc , ax ave , (c) T cool , ax ave and (d) D cool , ax ave .
Energies 15 05226 g029
Figure 30. Radial profile of P lin , rad ave (a) mean and (b) standard deviation.
Figure 30. Radial profile of P lin , rad ave (a) mean and (b) standard deviation.
Energies 15 05226 g030
Figure 31. Radial profile of T fc , rad ave (a) mean and (b) standard deviation.
Figure 31. Radial profile of T fc , rad ave (a) mean and (b) standard deviation.
Energies 15 05226 g031
Figure 32. ECDF of (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Figure 32. ECDF of (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Energies 15 05226 g032
Figure 33. Statistically significant cross-sections Spearman correlation coefficients for (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Figure 33. Statistically significant cross-sections Spearman correlation coefficients for (a) ρ i n , (b) P lin max , (c) T fc max and (d) H f max .
Energies 15 05226 g033
Figure 34. Temporal evolution of the (a) P lin , t max , (b) T fc , t max and (c) H f , t max .
Figure 34. Temporal evolution of the (a) P lin , t max , (b) T fc , t max and (c) H f , t max .
Energies 15 05226 g034
Figure 35. Radial profile of P lin , end ave (a) mean and (b) standard deviation.
Figure 35. Radial profile of P lin , end ave (a) mean and (b) standard deviation.
Energies 15 05226 g035
Figure 36. Radial profile of T fc , end ave (a) mean and (b) standard deviation.
Figure 36. Radial profile of T fc , end ave (a) mean and (b) standard deviation.
Energies 15 05226 g036
Figure 37. ECDF of (a) Bu ave and (b) Bu max .
Figure 37. ECDF of (a) Bu ave and (b) Bu max .
Energies 15 05226 g037
Figure 38. Statistically significant cross-sections Spearman correlation coefficients for the Bu max .
Figure 38. Statistically significant cross-sections Spearman correlation coefficients for the Bu max .
Energies 15 05226 g038
Figure 39. Evolution over depletion of the (a) Bu t max (b) P lin , t max and (c) T fc , t max .
Figure 39. Evolution over depletion of the (a) Bu t max (b) P lin , t max and (c) T fc , t max .
Energies 15 05226 g039
Figure 40. Axial profile of the (a) P lin , BOC ave (b) P lin , EOC ave and (c) Bu ax , EOC ave compared with the benchmark reference.
Figure 40. Axial profile of the (a) P lin , BOC ave (b) P lin , EOC ave and (c) Bu ax , EOC ave compared with the benchmark reference.
Energies 15 05226 g040
Figure 41. Radial profile of Bu rad , EOC ave (a) mean and (b) standard deviation at the EOC.
Figure 41. Radial profile of Bu rad , EOC ave (a) mean and (b) standard deviation at the EOC.
Energies 15 05226 g041
Figure 42. Average radial burnup discrepancy between the predicted P lin , BOC ave and the benchmark reference.
Figure 42. Average radial burnup discrepancy between the predicted P lin , BOC ave and the benchmark reference.
Energies 15 05226 g042
Table 1. Nuclide reactions analyzed in the uncertainty quantification framework.
Table 1. Nuclide reactions analyzed in the uncertainty quantification framework.
MTReaction
2Elastic scattering
4Inelastic scattering
18Fission
102Radiative capture
455 ν ¯ d : delayed neutrons released per fission
456 ν ¯ p : prompt neutrons released per fission
1018Fission spectrum
Table 2. CTF input uncertainty quantification and boundary conditions with N = Normal (Mean, Standard Deviation) and U = Uniform (Min, Max) determined based on LWR-UAM Phase II [37] and Phase III [24] specifications and recommendations.
Table 2. CTF input uncertainty quantification and boundary conditions with N = Normal (Mean, Standard Deviation) and U = Uniform (Min, Max) determined based on LWR-UAM Phase II [37] and Phase III [24] specifications and recommendations.
InputVariableProbability Density Function
Fuel Thermal Conductivity k f N (1.0, 0.00774)
Cladding Thermal Conductivity k c N (1.0, 0.00559)
Fuel Specific Heat Capacity C f N (1.0, 0.015)
Cladding Specific Heat Capacity C c N (1.0, 0.015)
Gap Conductance H gap U (0.75, 1.25)
Outlet Pressure p out N (1.0, 0.005)
Inlet Mass Flux m flux N (1.0, 0.001)
Inlet Coolant Temperature t in N (1.0, 0.00267)
Table 3. Manufacturing parameters input uncertainty quantification with N = Normal (Mean, Standard Deviation) determined based on LWR-UAM Phase II [37] and Phase III [24] specifications and recommendations.
Table 3. Manufacturing parameters input uncertainty quantification with N = Normal (Mean, Standard Deviation) determined based on LWR-UAM Phase II [37] and Phase III [24] specifications and recommendations.
InputVariableUnitsProbability Density Function
235U enrichment e uo 2 w/o%N (4.95, 0.000746)
UO 2 density ρ u o 2 g · cm 3 N (10.283, 0.0566)
UO 2 / Gd 2 O 3 density ρ g d g · cm 3 N (10.14, 0.0566)
Fuel outer radius rf o cmN (0.46955, 0.000433)
Cladding inner radius rc i cmN (0.4791, 0.0008)
Cladding outer radius rc o cmN (0.5464, 0.00083)
Instrumentation tube outer radius rit o cmN (0.6261, 0.00083)
Control rod outer radius rcr o cmN (0.6731, 0.00083)
Guide tube outer radius rgt o cmN (0.6731, 0.00083)
Table 4. Parameter coverage used for the history and branch modeling in the fuel assemblies.
Table 4. Parameter coverage used for the history and branch modeling in the fuel assemblies.
ParameterVariableUnitsMinMax
Fuel temperature T F K 5601600
Moderator density ρ M g · cm 3 0.608110.76971
Boron concentration C B ppm 01800
Control rod insertion--01
Table 5. Output quantities for steady-state calculations.
Table 5. Output quantities for steady-state calculations.
Output QuantityDescriptionVariable
ScalarMaximum linear power in the 3D core P lin max
Maximum fuel centerline temperature in the 3D core T fc max
Minimum coolant density in the 3D core D cool min
Effective multiplication factor k eff
Functional 1DRadially averaged axial linear power P lin , ax ave
Radially averaged axial fuel centerline temperature T fc , ax ave
Radially averaged axial coolant temperature T cool , ax ave
Radially averaged axial coolant density D cool , ax ave
Functional 2DAxially averaged radial linear power P lin , rad ave
Axially averaged radial fuel centerline temperature T fc , rad ave
Table 6. Output quantities for transient calculations.
Table 6. Output quantities for transient calculations.
Output QuantityDescriptionVariable
ScalarMaximum linear power during the transient P lin max
Maximum fuel centerline temperature during the transient T fc max
Maximum fuel enthalpy during the transient H f max
Inserted reactivity ρ i n
Functional 1DMaximum linear power over time P lin , t max
Maximum fuel centerline temperature over time T fc , t max
Maximum fuel enthalpy over time H f , t max
Functional 2DAxially averaged radial linear power, end of the transient P lin , end ave
Axially averaged radial fuel centerline temperature, end of the transient T fc , end ave
Table 7. Output quantities for depletion calculations.
Table 7. Output quantities for depletion calculations.
Output QuantityDescriptionVariable
ScalarAverage burnup at EOC Bu ave
Maximum burnup at EOC Bu max
1D FunctionalMaximum linear power over depletion P lin , t max
Maximum fuel centerline temperature over depletion T fc , t max
Maximum burnup over depletion Bu t max
Radially averaged axial linear power at BOC P lin , BOC ave
Radially averaged axial linear power at EOC P lin , BOC ave
Radially averaged axial burnup at EOC Bu ax , EOC ave
2D FunctionalAxially averaged radial burnup at EOC Bu rad , EOC ave
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Delipei, G.K.; Rouxelin, P.; Abarca, A.; Hou, J.; Avramova, M.; Ivanov, K. CTF-PARCS Core Multi-Physics Computational Framework for Efficient LWR Steady-State, Depletion and Transient Uncertainty Quantification. Energies 2022, 15, 5226. https://doi.org/10.3390/en15145226

AMA Style

Delipei GK, Rouxelin P, Abarca A, Hou J, Avramova M, Ivanov K. CTF-PARCS Core Multi-Physics Computational Framework for Efficient LWR Steady-State, Depletion and Transient Uncertainty Quantification. Energies. 2022; 15(14):5226. https://doi.org/10.3390/en15145226

Chicago/Turabian Style

Delipei, Gregory K., Pascal Rouxelin, Agustin Abarca, Jason Hou, Maria Avramova, and Kostadin Ivanov. 2022. "CTF-PARCS Core Multi-Physics Computational Framework for Efficient LWR Steady-State, Depletion and Transient Uncertainty Quantification" Energies 15, no. 14: 5226. https://doi.org/10.3390/en15145226

APA Style

Delipei, G. K., Rouxelin, P., Abarca, A., Hou, J., Avramova, M., & Ivanov, K. (2022). CTF-PARCS Core Multi-Physics Computational Framework for Efficient LWR Steady-State, Depletion and Transient Uncertainty Quantification. Energies, 15(14), 5226. https://doi.org/10.3390/en15145226

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop