Perturbation-Theory-Based Sensitivity and Uncertainty Analysis with CASMO-4

Abstract

The topic of this paper is the development of sensitivity and uncertainty analysis capability to the reactor physics code CASMO-4 in the context of the UAM (Uncertainty Analysis in Best-Estimate Modelling for Design, Operation, and Safety Analysis of LWRs) benchmark. The sensitivity analysis implementation is based on generalized perturbation theory, which enables computing the sensitivity profiles of reaction rate ratios efficiently by solving one generalized adjoint system for each response. Both the theoretical background and the practical guidelines for modifying a deterministic transport code to compute the generalized adjoint solutions and sensitivity coefficients are reviewed. The implementation to CASMO-4 is described in detail. The developed uncertainty analysis methodology is deterministic, meaning that the uncertainties are computed based on the sensitivity profiles and covariance matrices for the uncertain nuclear data parameters. The main conclusions related to the approach used for creating a covariance library compatible with the cross-section libraries of CASMO-4 are presented. Numerical results are given for a lattice physics test problem representing a BWR, and the results are compared to the TSUNAMI-2D sequence in SCALE 6.1.

1. Introduction

The topic of this paper is the development of sensitivity and uncertainty analysis capability to the reactor physics code CASMO-4 [1] in the context of the UAM (Uncertainty Analysis in Best-Estimate Modelling for Design, Operation and Safety Analysis of LWRs) benchmark [2]. At VTT, CASMO-4 is the standard tool for lattice physics calculations, and therefore it was a natural choice as the development platform for a sensitivity and uncertainty calculation system for the pin cell and fuel assembly exercises in the benchmark.

Sensitivities with respect to uncertain parameters can be computed efficiently by utilizing the adjoint system of the criticality equation. The propagated parameter uncertainty can then be calculated deterministically by the Sandwich rule by combining the sensitivity profiles with the covariance matrices of the parameters. As a first step, classical perturbation theory (CPT) was implemented to CASMO-4 to enable the computation of critical eigenvalue sensitivities with respect to nuclear data parameters. In this context, a methodology was devised for processing the covariance matrices from SCALE 6 [3] to become compatible with the cross-section libraries of CASMO-4 to enable uncertainty analysis. This work has been reported in detail in [4]. Recently, generalized perturbation theory (GPT) has been added to the code as a new feature. This enables performing sensitivity analysis for responses that can be presented as reaction rate ratios. In this framework, one generalized adjoint system needs to be solved for each response, after which the response sensitivity profiles for all parameters of interest can be computed in an efficient manner.

This paper is organized as follows. Section 2 reviews the theoretical background for sensitivity and uncertainty analysis based on generalized perturbation theory, and Section 3 focuses on the implementation to CASMO-4. In Section 3.1, the computation of generalized adjoint solutions is considered and practical guidelines are presented for modifying a deterministic transport code to solve the adjoint problems needed in sensitivity analysis. Section 3.2 concerns the computation of sensitivity and uncertainty profiles. Finally, in Section 4, numerical results are presented for a lattice physics test problem representing a BWR, and they are compared to the TSUNAMI-2D sequence in SCALE 6.1.

2. Theoretical Background

The purpose of sensitivity analysis is to study how sensitive a mathematical model is to perturbations in its uncertain parameters. The target of uncertainty analysis is to estimate how the uncertainty in these parameters is propagated to a response dependent on the mathematical model under consideration. In this work the mathematical model is the neutron transport eigenvalue problem, which can be written in operator form as 1𝐀Φ=𝑘𝐁Φ, (1) where Φ∈𝐻_Φ is the neutron flux, 𝐻_Φ is a Hilbert space, and 𝑘 is the multiplication factor. The uncertain parameters consist of nuclear data parameters and they are denoted by the vector 𝝈∈𝐸_𝜎. It should be noted that both the continuous-energy criticality equation and the various systems derived from it in numerical computations can be written in the form of (1).

2.1. Sensitivity Analysis

The object of local sensitivity analysis is to determine how the response 𝑅 depends on the uncertain parameters in the vicinity of their best-estimate values. In this work, the responses under consideration include homogenized assembly parameters and the multiplication factor, whereas the uncertain parameters are neutron cross-sections. When considering the continuous-energy eigenvalue problem, the cross-sections are functions of energy and location, and the appropriate derivative is the functional directional derivative called the Gâteaux variation [5]. It follows that the sensitivity of 𝑅 with respect to the perturbation 𝐡=[𝛿Φ,𝛿𝝈]∈𝐷=𝐻_Φ×𝐸_𝜎 at the point 𝐞̂=[Φ,𝝈]∈𝐷 may be defined as 𝛿𝑅(𝐞̂;𝐡)=lim_𝑡→0𝑅(𝐞̂+𝑡𝐡)−𝑅(𝐞̂)𝑡. (2)

When the parameters 𝝈 are perturbed, also the solution Φ is affected and therefore the computation of the sensitivity 𝛿𝑅(𝐞̂;𝐡) requires that the perturbation 𝛿Φ is known. In principle, 𝛿Φ can be computed to first order from the following forward sensitivity system: (̂)1𝛿𝐀𝐞;𝐡=−𝑘²(̂)1𝛿𝑘𝐞;𝐡𝐁Φ+𝑘𝛿𝐁(𝐞̂;𝐡)⟺𝐀^′_𝜎(𝐞̂)(𝐞̂)1𝛿𝝈+𝐀𝛿Φ=−𝑘²(̂)+1𝛿𝑘𝐞;𝐡𝐁Φ𝑘𝐁^′_𝜎(𝐞̂)1𝛿𝝈+𝑘𝐁(𝐞̂)𝛿Φ, (3) which can be derived by taking the Gâteaux variation of system (1) with respect to a perturbation 𝐡 on both sides. However, when computing several sensitivities, this approach would require the repetitive solving of (3).

Fortunately, the sensitivities can be computed more efficiently by exploiting the adjoint of (1), which is defined as the system that satisfies the following relation: (In some cases the adjoint relation needs to be written in the form ⟨𝐀Φ+(1/𝑘)𝐁Φ,Ψ⟩=⟨Φ,𝐀^∗Ψ+(1/𝑘)𝐁^∗Ψ⟩+[𝐏(Ψ,Φ)]_{𝐱∈𝜕Ω}, where [𝐏(Ψ,Φ)]_{𝐱∈𝜕Ω} is a bilinear form associated with the system. We will only consider cases where it is straightforward to force this term to vanish.) 1𝐀Φ−𝑘𝐁Φ,Ψ=Φ,𝐀^∗1Ψ−𝑘𝐁^∗Ψ, (4) where the brackets ⟨⋅,⋅⟩ denote an inner product. When considering the continuous-energy criticality equation, it is customary to employ the 𝐿² inner product [6, 7]. The solution to the adjoint problem 𝐀^∗−1𝑘𝐁^∗Ψ=0 (5) is called the fundamental adjoint. Physically, the solution to this system can be interpreted to represent the average contribution, that is, importance of a neutron to the multiplication factor. Interestingly, the adjoint system of (5) can be derived solely based on this physical interpretation [8]. Like the neutron flux, the fundamental adjoint has an arbitrary normalization, and the concept of importance should be understood in relative terms. Therefore, the value Ψ(𝐫,Ω,𝐸) represents the importance of a neutron located at the point [𝐫,Ω,𝐸] compared to the importance of neutrons elsewhere in the phase space [9]. Based on this physical reasoning, it can be deduced that the fundamental adjoint must always be nonnegative.

By utilizing (4) and (5), it is straightforward to obtain the following expression for the relative sensitivity of the multiplication factor with respect to a perturbation 𝛿𝝈 (For derivation, see for example, [4, 10]): 𝛿𝑘(𝐞̂;𝐡)𝑘𝐀=−^′𝜎(𝐞̂)−(1/𝑘)𝐁^′_𝜎(𝐞̂)𝛿𝝈,Ψ)⟩.⟨(1/𝑘𝐁Φ,Ψ (6) This equation is known in reactor physics as classical perturbation theory. In addition, the adjoint system can be utilized in the sensitivity analysis of the eigenvalue problem for other responses fulfilling the following properties. Firstly, the response 𝑅 must be Fréchet-differentiable with respect to Φ, in which case we can write 𝛿𝑅(𝐞̂;𝐡)=𝑅^′_𝜎(𝐞̂)𝛿𝝈+∇_Φ𝑅(𝐞̂),𝛿Φ_Φ. (7) In addition, the response's Fréchet derivative ∇_Φ𝑅 (also called gradient) must be orthogonal to the forward solution ∇_Φ𝑅,Φ=0. (8) When these assumptions are fulfilled, the generalized adjoint corresponding to the response 𝑅 can be defined as the solution to the following inhomogeneous system: 𝐀^∗+1𝑘𝐁^∗∇Γ=^Φ𝑅𝑅. (9) Notice that in the previous equation the eigenvalue 𝑘 is fixed to correspond to the solution of (1) and therefore the operator 𝐀^∗+(1/𝑘)𝐁^∗ is singular, which necessitates (8) in order for the solution Γ to exist. Also, when a solution Γ₀ to (9) exists, there exists an infinite amount of solutions of the form Γ=Γ₀+𝑎Ψ,𝑎∈ℝ. (10) In this case, it is possible to choose a solution that is orthogonal to the (forward) fission source. This particular solution can be written as Γ_𝑝=Γ₀−⟨Γ₀,𝐁Φ⟩Ψ⟨Ψ,𝐁Φ⟩=Γ₀−⟨𝐁^∗Γ_𝟎,Φ⟩⟨𝐁^∗Ψ,Φ⟩Ψ. (11)

We can now derive a practical expression for the response sensitivity with respect to a perturbation 𝛿𝝈: 𝛿𝑅(𝐞̂,𝐡)𝑅=𝑅^′_𝜎(𝐞̂)𝛿𝝈𝑅+∇_Φ𝑅(𝐞̂)𝑅,𝛿Φ_Φ(9)=𝑅^′_𝜎(𝐞̂)𝛿𝝈𝑅+(𝐀^∗+1𝑘𝐁^∗)Γ,𝛿Φ_Φ(4)=𝑅^′_𝜎(𝐞̂)𝛿𝝈𝑅+1Γ,𝐀+𝑘𝐁𝛿Φ_Φ(3)=𝑅^′_𝜎(𝐞̂)𝛿𝝈𝑅−Γ,𝐀^′_𝜎(𝐞̂)−1𝑘𝐁^′_𝜎(𝐞̂)𝛿𝝈_Φ−𝛿𝑘(𝐞̂;𝐡)𝑘²⟨Γ,𝐁Φ⟩_Φ=𝑅^′_𝜎(𝐞̂)𝛿𝝈𝑅−Γ_𝑝,𝐀^′_𝜎(𝐞̂)−1𝑘𝐁^′_𝜎(𝐞̂)𝛿𝝈_Φ. (12) This framework is often referred to as generalized perturbation theory when the response 𝑅 is of the form: 𝑅(𝐞)=Φ,𝚺₁Φ,𝚺₂. (13) In this case, it is straightforward to show that (8) is satisfied and that 𝑅 is Fréchet-differentiable, the relative gradient being ∇_Φ𝑅𝑅=𝚺₁Φ,𝚺₁−𝚺₂Φ,𝚺₂. (14)

The generalized adjoint Γ(𝐫,Ω,𝐸) can be physically interpreted as the average contribution of an additional neutron at the phase space point [𝐫,Ω,𝐸] to the response under consideration. The generalized adjoint is normalized according to the value of the response. It should also be noticed that since an additional neutron may also reduce the value of the response, generalized adjoints can also have negative values. The gradient of the response may also be negative in some parts of the phase space.

In practice, the eigenvalue problem and the corresponding adjoint equations are solved numerically, which gives rise to some complications in the perturbation theory formalism. Ideally, the discretizations should be performed in a consistent manner so that the respective adjoint relations are satisfied at all stages of the computation [5]. However, this is usually infeasible in reactor physics calculations, and therefore it is customary to take the eigenvalue problem discretized with respect to energy and direction as the starting point for sensitivity analysis. This issue is discussed in more detail in [4].

2.2. Uncertainty Analysis

The uncertainty of the uncertain parameters 𝝈 should be understood in terms of the Bayesian probability interpretation. In this framework, probability is defined as a subjective measure that characterizes the plausibility of various hypotheses. When estimating parameters, all knowledge about a parameter 𝜎_𝑗 is assumed to be incorporated into its marginal probability distribution 𝑝(𝜎_𝑗). This distribution is defined so that the integral ∫^𝑏_𝑎𝑝(𝜎_𝑗)𝑑𝜎_𝑗 corresponds to the (Bayesian) probability that the value of 𝜎_𝑗 belongs to the interval [𝑎,𝑏]. The distribution 𝑝(𝝈) can then be used to form an estimate for the parameters and their associated uncertainties. Usually, the mean value or the mode is chosen as the estimate for the parameters, whereas the covariance matrix of the distribution is chosen as the descriptive statistic for the uncertainty.

In the Bayesian formalism, the outcome of the uncertainty analysis should ideally be the full posterior distribution 𝑝(𝐑). However, determining 𝑝(𝐑) analytically is usually extremely challenging and the distribution can only be estimated pointwise based on a simulation. In deterministic uncertainty analysis, the objective is not to form the entire distribution 𝑝(𝐑), but to compute an estimate for the covariance matrix Cov[𝐑] by linearizing the response 𝐑≈𝐒𝝈. Here 𝐒∈ℝ^𝐽×𝐾 is the response vector sensitivity matrix, 𝐽 is the number of responses, and 𝐾 is the number of uncertain parameters. After linearizing the response, the covariance matrix can be computed simply using the identity Cov[𝐑]≈Cov[𝐒𝝈]=𝐒Cov[𝝈]𝐒^𝑇, (15) known as the first-order uncertainty propagation formula or the Sandwich rule.

3. Implementation

3.1. Computation of Generalized Adjoint Fluxes

This section reviews the guidelines for modifying a deterministic transport solver to compute the adjoint solutions needed in generalized perturbation theory and describes the methodology used in the implementation to CASMO-4. As mentioned previously, the description on the implementation of classical perturbation theory to CASMO-4 has been recently published in [4], and therefore, in this paper, the emphasis is placed on the GPT-specific features.

As explained in Section 2.1, it is customary to take the energy- and direction-discretized system as the starting point for perturbation theory. In CASMO-4, the multigroup criticality equation is solved by the method of characteristics assuming isotropic scattering. Therefore, the following system of equations may be taken as the forward problem: Ω_𝑚⋅∇Φ^𝑔𝐫,Ω_𝑚+Σ^𝑔Φ^𝑔𝐫,Ω_𝑚=14𝜋^𝐺_ℎ=1Σ^ℎ→𝑔_s𝜙^ℎ(𝐫)+𝜒_𝑔4𝜋𝑘^𝐺_ℎ=1𝜈Σ^ℎ_f𝜙^ℎ(𝐫),𝑔=1,…,𝐺. (16) In (16) the scalar flux is approximated by the quadrature formula 𝜙^ℎ(𝐫)=^𝑀_𝑚=1𝜔_𝑚Φ^ℎ𝐫,Ω_𝑚. (17) In order to simulate an infinite lattice, the boundary conditions are often assumed to be reflective, that is, Φ𝐫,Ω_𝑚,𝐸=Φ𝐫,Ω^′_𝑚,𝐸,𝐫∈Γ,Ω_𝑚⋅𝐧<0, (18) where Ω_𝑚=Ω^′_𝑚−2(𝐧⋅Ω^′_𝑚)𝐧 is the reflection direction. The inner product corresponding to this discretization can be defined in a consistent manner as ⟨Φ,Ψ⟩=^𝐺𝑀_𝑔=1𝑚=1𝜔_𝑚𝐷𝑑³𝐫Φ^𝑔𝐫,Ω_𝑚Ψ^𝑔𝐫,Ω_𝑚. (19) The adjoint system can now be written −Ω_𝑚⋅∇Ψ^𝑔𝐫,Ω_𝑚+Σ^𝑔Ψ^𝑔𝐫,Ω_𝑚=14𝜋^𝐺_ℎ=1Σ^𝑔→ℎ_s𝜓^ℎ(𝐫)+𝜈Σ^𝑔_f4𝜋𝑘^𝐺_ℎ=1𝜒_ℎ𝜓^ℎ(𝐫),𝑔=1,…,𝐺, (20) with the boundary conditions Ψ𝐫,Ω_𝑚,𝐸=Ψ𝐫,Ω^′_𝑚,𝐸,𝐫∈Γ,Ω_𝑚⋅𝐧>0. (21) It is straightforward to check that the systems (16) and (20) with their respective boundary conditions satisfy (4) with respect to the inner product defined by (19).

The generalized adjoint problem for a response of the form of (13) can now be written −Ω_𝑚⋅∇Γ^𝑔𝐫,Ω_𝑚+Σ^𝑔Γ^𝑔𝐫,Ω_𝑚=14𝜋^𝐺_ℎ=1Σ^𝑔→ℎ_s𝛾^ℎ(𝐫)+𝜈Σ^𝑔_f4𝜋𝑘^𝐺_ℎ=1𝜒_ℎ𝛾^ℎ(𝐫)+Σ^𝑔₁(𝐫)Φ,𝚺₁−Σ^𝑔₂(𝐫)Φ,𝚺₂,𝑔=1,…,𝐺, (22) where the generalized adjoint of the scalar flux has been denoted by 𝛾^ℎ(𝐫). As explained in Section 2.1, this system may have an infinite number of solutions, of which we wish to solve the one that satisfies 𝐁^∗Γ_𝑝,Φ=0. (23)

In deterministic transport solvers, the iteration for fixed source calculations is generally of the form 𝐀Φ^𝑛+1=𝐁Φ^𝑛+𝐒, (24) where 𝐒 is an external source. This iteration scheme with a fixed eigenvalue is also well suited for solving the generalized adjoint problem of (22), in which case the iteration takes the form 𝐀^∗Γ^𝑛+1=1𝑘𝐁^∗Γ^𝑛+∇_Φ𝑅𝑅. (25) During the iteration, however, the convergence to the particular solution that is orthogonal to the fission source must be ensured. It is straightforward to show that if the initial guess for the generalized adjoint flux satisfies (23), this orthogonality property is preserved during the iteration. Firstly, 𝐀^∗Γ^𝑛+1,Φ=⁽⁾²⁵1𝑘𝐁^∗Γ^𝑛+∇,Φ^Φ𝑅𝑅,Φ₍₈₎=1𝑘𝐁^∗Γ^𝑛,Φ. (26) On the other hand, 𝐀^∗Γ^𝑛+1,Φ=Γ^𝑛+1,𝐀Φ=⁽¹⁾1𝑘Γ^𝑛+1=1,𝐁Φ𝑘𝐁^∗Γ^𝑛+1,Φ. (27) Therefore, for each iteration 𝑛, 𝐁^∗Γ^𝑛+1,Φ=𝐁^∗Γ^𝑛,Φ, (28) from which the result follows. In practice, however, due to round-off errors and the unavoidable inconsistencies in formulating the discretizations and adjoint relations, a refinement of the iteration scheme is necessary to guarantee that (23) remains satisfied [11]. A suitable procedure is to force the orthogonality of the solution with each outer iteration. In this case, in accordance with (11), the iteration takes the form 𝐀^∗Γ^𝑛+1=1𝑘𝐁^∗Γ^𝑛−⟨𝐁^∗Γ^𝑛,Φ⟩⟨𝐁^∗Ψ+∇Ψ,Φ⟩^Φ𝑅𝑅. (29) Notice that this iteration scheme requires that the forward solution and the fundamental adjoint solution have been previously computed and that they are accessible during the iteration.

By comparing (29) with the forward problem of (16), it can be seen that if the forward system had an external source, the systems would be of the same form with the exception that the adjoint system is solved in the opposite direction. Therefore, if the transport solver does not rely on the assumption of the nonnegativity of the flux or the sources, relatively few modifications are needed to transform the solver to also compute the generalized adjoint functions. For example, the method of characteristics, used in CASMO-4, does not require that the solution or the sources are non-negative. In this case, the following operations need to be performed before the adjoint calculation [10].(1)Transpose the scattering matrix. (2)Interchange the vectors 𝝂𝝈_f and 𝝌. (3)Invert the group indices for all variables as follows: 𝐺↔1, (𝐺−1)↔2,….After these operations, the transport solver can be used to compute the fundamental adjoint solution. Notice also that these operations automatically convert the forward boundary conditions to the adjoint boundary conditions. When solving a generalized adjoint problem, the following changes need to be additionally implemented within the (forward) transport solver.(1)Add the response gradient ∇_Φ𝑅/𝑅 to the variable for an external source. (2)Modify the fission source 𝐹_𝑔 to the form 𝐹_𝑔=𝜒_𝑔4𝜋𝑘^𝐺_ℎ=1𝜈Σ^ℎ_f𝜙^ℎ(𝐫)−𝐁Φ,Φ_F𝐁Ψ,Φ_FΨ^𝑔, (30) where Φ_F denotes the forward solution of (16) and Ψ the adjoint solution of (20).

The multigroup solution Φ given by the solver must then be interpreted so that, for example, Φ^𝑔(𝐫,Ω) corresponds to Γ^{𝐺+1−𝑔}(𝐫,−Ω). Notice that if the transport solver is based on a numerical scheme that relies on the nonnegativity of the flux or the sources, some additional modifications are necessary in addition to the ones described above. For further details, see for example, [11].

3.2. Computation of Sensitivity and Uncertainty Profiles

After obtaining the adjoint solutions, the sensitivities with respect to the multigroup nuclear data parameters can be computed according to (6) and (12). Notice that even after the multigroup approximation, these parameters are still spatial functions and therefore the derivatives in the equations refer to functional derivatives. The inner product in the sensitivity expressions can be discretized as ⟨Φ,Ψ⟩≈^𝐼𝐺_𝑖=1^𝑀_𝑔=1𝑚=1𝜔_𝑚𝑉_𝑖Φ^{𝑔,𝑖,𝑚}Ψ^{𝑔,𝑖,𝑚}, (31) where 𝑖 denotes the mesh index and Φ^{𝑔,𝑖,𝑚} and Ψ^{𝑔,𝑖,𝑚} denote the average fluxes.

In order to compute the uncertainties using the Sandwich rule, the sensitivities and covariance matrices need to be formed with respect to the same parameters using the same energy group structure. In the SCALE 6 covariance library [3], the available covariance matrices are given in a 40-group structure for the parameters listed in Table 1. Most of these covariance matrices are nuclide specific. It should be emphasized that there is no covariance data for the group-to-group transfer cross-sections.

Table 1 

Parameters for which there exists covariance data in the SCALE library.

Multigroup covariance matrices can in principle be transformed to another multigroup structure by simple mathematical techniques. The applicability of this approach depends on the differences between the group structures. In particular, the widths of the energy groups should not dramatically change. In this work, the code Angelo 2.3 [12] was used to transform the matrices to the energy group structure used in the sensitivity calculations with CASMO-4. The transformation procedure used in the code is based on flat-flux approximation, where the resampled values on the new grid are computed as lethargy overlap weighted averages. For further details, see [13]. When modifying the energy group structure of fission spectrum covariance matrices, further correction procedures are necessary in order to guarantee that the covariance matrices are in accordance with the normalization condition ∑_𝑔𝜒^𝑔=1 [14]. The correction can also be applied to the fission spectrum sensitivities in which case the sensitivities are called constrained [14]. This was the approach chosen in this work.

In order to utilize the covariance data given for the parameters in Table 1, sensitivity profiles should be computed with respect to the same parameters. However, many lattice physics codes such as CASMO, HELIOS [15], WIMS [16], and DRAGON [17] employ nuclear data libraries that do not contain cross-section data for the individual capture and scattering reactions, but only for the total capture and scattering cross-section. There are generally three different approaches to overcome this difficulty. The most natural approach is perhaps to add the missing cross-sections to the code, either by creating a new cross-section library or by modifying the cross-sections inside the code [18]. Another option, suitable for deterministic analysis, is not to use problem-dependent cross-sections in the sensitivity analysis. In this case, the sensitivity coefficients can be computed outside the code based on the forward and adjoint fluxes and any set of cross-sections. This was the idea, for example, behind connecting DRAGON with the sensitivity and uncertainty analysis code SUSD3D after a generalized adjoint mode was implemented to DRAGON [19]. The third option is to form the covariance matrices corresponding to the total capture and scattering cross-sections [4]. This is the approach that was chosen in this work.

Since the relationships between the total and individual capture and scattering reactions are linear, the covariance matrices corresponding to the total capture and scattering reactions can be computed with the Sandwich rule without introducing any approximation. The method used for combining the covariance matrices has been recently described in detail in [4]. Therefore, only the most important conclusions related to the methodology are repeated here.

Firstly, in the context of the capture reactions, the results are expected to be fully consistent with the case where the sensitivities are computed with respect to the individual capture reactions. In the case of the scattering reactions, however, the sensitivity profiles with respect to the individual and the total scattering cross-sections cannot be defined in a consistent manner and this affects the uncertainty results. In this context, it should be emphasized that the treatment of the covariance matrices involves no approximations and the inconsistency is solely related to the computation of the sensitivities. As mentioned previously, there is no cross-section data for the transfer cross-sections 𝜎^{ℎ→𝑔,𝑗}_x but only for 𝜎^𝑔,𝑗_x=∑^𝐺_ℎ=1𝜎^{𝑔→ℎ,𝑗}_x, where x refers to a scattering reaction (e.g., elastic, inelastic) and 𝑗 is the nuclide index. Therefore, in order to use the scattering covariance data, the sensitivity profiles should be computed with respect to 𝜎^𝑔,𝑗_x. Because of the scattering source term in (16), however, the derivative with respect to 𝜎^𝑔,𝑗_x is not mathematically well defined without additional constraints. Typically it is assumed that the probabilities of transfers to various groups are fixed, that is, 𝜎^{𝑔→ℎ,𝑗}_x=𝜎^𝑔,𝑗_x𝑝^{𝑔→ℎ,𝑗}_x, (32) where 𝑝^𝑔→ℎ_x is the proportion of neutrons scattered from energy group 𝑔 to energy group ℎ, which is assumed to remain fixed even if the scattering cross-section 𝜎^𝑔,𝑗_x is perturbed [20]. Based on this assumption, the scattering source in (16) can be written as 𝑆^𝑔=14𝜋^𝐺_ℎ=1Σ^ℎ→𝑔_s𝜙^ℎ=14𝜋_x𝑗𝑁^𝑗𝐺_ℎ=1𝜎^ℎ,𝑗_x𝑝^ℎ→𝑔_x𝜙^ℎ, (33) where the summations over x include all scattering reactions. After this assumption, the derivative with respect to 𝜎^𝑔,𝑗_x is well defined and can be computed as usual. It is straightforward to show that this approach corresponds to computing the sensitivity coefficients with respect to the transfer cross-sections 𝜎^{𝑔→ℎ,𝑗}_x and summing them over ℎ.

However, the sensitivity with respect to the total scattering cross-section 𝜎^𝑗_s=∑_𝑥𝜎^𝑗_x is not well defined if the constraint (32) is enforced. In order to define this sensitivity, fixed transfer rates must be assumed for the total scattering cross-section. Also, computing the total scattering sensitivity as the sum of the individual scattering sensitivities implicitly enforces this constraint. Since the two assumptions required to compute the individual and total scattering sensitivities are inconsistent, the chain rule of derivation does not apply to them, and, for example, although 𝜎^𝑔,𝑗_s=𝜎^𝑔,𝑗_e+𝜎^𝑔,𝑗_i holds, 𝑑𝑅/𝑑𝜎^𝑔,𝑗_e≠(𝑑𝑅/𝑑𝜎^𝑔,𝑗_s)(𝑑𝜎^𝑔,𝑗_s/𝑑𝜎^𝑔,𝑗_e).

4. Numerical Results for PB-2 Lattice Physics Exercise

The calculation framework was applied to the BWR test case from the UAM benchmark lattice physics Exercise 1.2 considering a single fuel assembly with reflective boundary conditions [2], and the results were compared against the TSUNAMI-2D sequence in SCALE 6.1 [21]. The test problem represents Peach Bottom 2 (PB-2) under hot zero power conditions. Two-group homogenized cross-sections have been considered as responses in the GPT framework.

The outline of the CASMO-4 calculations is presented in Figure 1. The calculations were carried out using the cross-section library E60200 that contains 70 energy groups and is based on ENDF/B-VI data [22]. The covariance data were taken from the SCALE 6 library ZZ-SCALE6.0/COVA-44G [3] according to the guidelines of the benchmark. The library is based on evaluations from various sources (including ENDF/B-VII, ENDF/B-VI, JENDL-3.1) and approximate covariance data. The covariances in the library are given in relative terms, and therefore the library is intended to be used with all cross-section libraries including the ones that are inconsistent with the evaluations. While this is not strictly correct, it is considered to be acceptable due to the scarcity of comprehensive covariance data among other reasons [23].

Figure 1 

Outline of the CASMO-4 calculations.

The list of the nuclides present in these test cases can be found in the benchmark specification [2]. Apart from the isotopes of chromium and iron, all available covariance data in the library was included in the uncertainty computations. The reason for excluding these isotopes is that the employed cross-section library E60200 does not contain isotope-specific cross-sections for these materials but only cross-sections for natural chromium and iron.

The covariance matrices from ZZ-SCALE6.0/COVA-44G were processed for compatibility with CASMO-4. The sensitivity profiles in CASMO-4 were computed using the 40-group structure option that was the closest match to the amount of groups in the covariance data and, as mentioned in Section 3.2, the code Angelo 2.3 [12] was used to process the covariance matrices to this energy group structure. Next, the nuclear data processing code NJOY [24] was used to transform the 40-group covariance files to the BOXR format. Auxiliary FORTRAN programs were written for combining the covariance matrices according to the principles described in Section 3.2.

The TSUNAMI-2D calculations were performed using the ENDF/B-VI-based cross-section library V6-238 containing 238 energy groups. The module CENTRM was used for self-shielding. Implicit sensitivity analysis [9] was omitted in the TSUNAMI calculations in order to facilitate the comparison of the results given by CASMO-4 and TSUNAMI-2D.

4.1. Results Based on Classical Perturbation Theory

A summary of the results based on classical perturbation theory for the multiplication factor is presented in Table 2. The relative difference between the multiplication factors computed with CASMO-4 and TSUNAMI-2D is 52 pcm in both forward and adjoint cases. For the total uncertainty, the values given by CASMO-4 and TSUNAMI-2D are also very consistent. Table 3 shows the five most significant sources of uncertainty together with the corresponding energy- and region-integrated sensitivity coefficients. As can be seen from this table, both the sensitivity and the uncertainty results are in good accordance. The greatest difference occurs for the capture cross-section of ²³⁸U, for which CASMO-4 yields a greater sensitivity. This appears to originate from the differences in the cross-section libraries. In particular, the cross-section library E60200 used in the CASMO-4 calculation has not been reduced in terms of the ²³⁸U resonance integral, which is known to be overestimated in the ENDF/B-VI data [22].

Table 2 

Summary of the results for the multiplication factor.

Table 3 

The five most significant sources of uncertainty for the multiplication factor and the corresponding energy- and region-integrated relative sensitivity coefficients.

Figure 2 shows the volume-averaged forward flux and the volume-averaged fundamental adjoint Ψ corresponding to this test case. As explained in Section 2.1, the value Ψ^𝑔 represents the average importance of neutrons in the energy group 𝑔 to the multiplication factor in comparison to neutrons in other energy groups. The plot can be interpreted from this point of view. For example, it can easily be seen from the figure how the ²³⁸U capture cross-section resonances reduce the importance of neutrons in the corresponding energy groups. This phenomenon is particularly clear in the energy group 𝐸₁₄=[4.00eV,9.88eV], where the multigroup capture cross-section attains its maximum value. It can also be clearly distinguished from the plot how the adjoint function has a higher value in the energy groups corresponding to the peaks in the fission cross-section of ²³⁵U. To further demonstrate this, Figure 3 shows a plot of the problem-dependent ²³⁵U fission and ²³⁸U capture cross-sections in the same 40-group structure. The increase in the adjoint values in the highest energy groups corresponds mainly to the increase in the value of 𝜈 at these energies.

Figure 2 

Volume-averaged forward flux and fundamental adjoint flux.

Figure 3 

Self-shielded multigroup cross-sections corresponding to the test problem.

Figure 4 shows the multiplication factor sensitivity profiles for the parameters, whose integrated sensitivity coefficients have the greatest absolute values, excluding the sensitivity profile with respect to the fission spectrum of ²³⁵U, which was constrained in the computation. As can be seen from the figure, the multiplication factor is the most sensitive to the fission parameters of ²³⁵U, the capture cross-section of ²³⁸U, and the scattering cross-section of ¹H. The positive sensitivity to the capture of ²³⁸U in the highest energy group follows from the fact that in CASMO-4 the (n,2n) reaction cross-section has been included in the capture cross-section with a negative sign in this group. It is instructive to compare the sensitivity profiles with the forward and adjoint fluxes plotted in Figure 2. Notice that the peaks in the sensitivity profiles of ²³⁵U coincide with the thermal peak of the neutron flux, where most of the fissions occur. In general, perturbing a nuclear parameter has a greater impact on the results in the energy groups, where the flux is higher. On the contrary, the values of the fundamental adjoint represent the average importance of neutrons in different energy groups. In particular, the lowest energy group has the highest importance, but this is not manifested in the sensitivity profiles, as the flux is very close to zero in this group. The negative sensitivities to the scattering reaction of ¹H in the four lowest energy groups can be attributed to the fact that in these groups upscattering is more likely than downscattering. Therefore, neutrons are scattered to energy groups with a lower importance. The same reasoning applies to the scattering sensitivity of ¹H in the highest energy groups, where neutrons are scattered downwards and the values of the adjoint function decrease rapidly with energy.

Figure 4 

Multiplication factor sensitivity profiles.

4.2. Results Based on Generalized Perturbation Theory

Table 4 presents the values and the total uncertainties of the homogenized two-group cross-sections that were considered as responses in the GPT-based sensitivity and uncertainty analysis. In computing the responses, the thermal cut-off was set at 0.625 eV. It can be seen from the table that all total uncertainty values are in good agreement with the thermal responses, whereas for the fast responses the uncertainties given by TSUNAMI-2D are consistently greater.

Table 4 

Values and uncertainties of the responses considered in the GPT framework.

Tables 5 and 6 show more detailed sensitivity and uncertainty results for the two-group homogenized production cross-sections 𝜈Σ_f,1 and 𝜈Σ_f,2. As can be seen from Table 6, in the case of 𝜈Σ_f,1, the difference in the total uncertainty values given by CASMO-4 and TSUNAMI-2D is attributable to the scattering of ²³⁸U, for which TSUNAMI-2D yields a significantly greater uncertainty value, although the total scattering sensitivity coefficients given by both codes are very close. As explained in Section 3.2, the sensitivity with respect to the total scattering cross-section can only be defined if the group-to-group transfer probabilities are assumed to be fixed for the total scattering. Also, defining the total scattering sensitivity as the sum of the individual scattering sensitivities implicitly enforces this assumption. However, in the TSUNAMI-2D computation, the total scattering uncertainty is computed based on the individual scattering sensitivities, which rely on the assumption of fixed transfer rates for each scattering reaction. The difference in the total scattering uncertainties is hence explained by incompatible constraints in the two uncertainty calculations. This phenomenon is more evident for the fast group responses since they are more sensitive to the inelastic scattering of ²³⁸U.

Table 5 

The five most significant sources of uncertainty for the response 𝜈Σ_f,2 and the corresponding energy- and region-integrated relative sensitivity coefficients.

Table 6 

The five most significant sources of uncertainty for the response 𝜈Σ_f,1 and the corresponding energy- and region-integrated relative sensitivity coefficients. The sensitivity coefficients with respect to the parameter 𝜒 have been constrained.

Figure 5 shows the volume-averaged generalized adjoint solutions for the responses 𝜈Σ_f,1 and 𝜈Σ_f,2, denoted by Γ_f,1 and Γ_f,2, respectively. As previously explained, the adjoint values in each energy group can be interpreted to represent the average importance of neutrons in that group to the response under consideration. Therefore, it is not surprising that thermal neutrons are more important to the response 𝜈Σ_f,2, whereas fast neutrons are more important to the response 𝜈Σ_f,1. The positive values of Γ_f,2 in the fast groups result from the downscattering of neutrons. Notice that Γ_f,1 also has a small positive value in the first thermal group, which corresponds to the possibility of upscattering. For the most part, both adjoint fluxes qualitatively follow the fission cross-section of  ²³⁵U plotted in Figure 3. In the highest energy groups, the values of Γ_f,1 increase rapidly due to the increase in the values of 𝜈. The negative values of Γ_f,1 between 0.111 MeV and 2.231 MeV signify that additional neutrons in those energy groups would on average contribute more to the denominator ⟨Φ,𝟏⟩₁ than to the numerator ⟨Φ,𝝂𝚺_f⟩₁. This in accordance with the fact that fission is unlikely to occur in this energy region.

Figure 5 

Volume-averaged generalized adjoints corresponding to the responses 𝜈Σf,1 and 𝜈Σf,2.

Figure 6 shows the sensitivity profiles of 𝜈Σ_f,1 with respect to the parameters, whose integrated sensitivity coefficients have the greatest absolute values. As can be anticipated, the response is the most sensitive to the fission parameters of ²³⁵U and ²³⁸U and in addition to the scattering of ¹H. It is interesting to compare these profiles with the plot of the generalized adjoint Γ_f,1 in Figure 5. The sensitivity to the scattering of ¹H has the smallest values in the groups with the highest importance, as this reaction transfers neutrons to energy groups with a lower importance. Since fast neutrons mostly scatter downwards, the scattering sensitivity has positive values in the groups between 149 eV and 1.35 MeV, where the importance decreases with increasing energy. This trend is reversed at 1.35 MeV, where the importance of the energy groups begins to increase with energy, mainly due to the increase in the values of 𝜈 at these energies.

Figure 6 

Sensitivity profiles of the response 𝜈Σf,1.

The sensitivity profiles corresponding to the response 𝜈Σ_f,2 are plotted in Figure 7. It is noteworthy that the profiles qualitatively resemble the respective profiles of the multiplication factor in the thermal region, whereas they quickly fall to nearly zero in the fast region. From the perspective of the GPT framework, it is again enlightening to compare the sensitivity plots with the adjoint function Γ_f,2 plotted in Figure 5 and the flux Φ shown in Figure 2. In the case of this response, the average importance of neutrons increases steadily with decreasing energy. Therefore, it is reasonable that the scattering sensitivities are again negative in the groups where upscattering is more likely than downscattering. Also, the sensitivities peak in the energy region coinciding with the thermal peak of the forward flux.

Figure 7 

Sensitivity profiles of the response 𝜈Σf,2.

5. Summary and Conclusions

Sensitivity and uncertainty analysis capability has been developed to the reactor physics code CASMO-4 in the context of the UAM benchmark. Sensitivities with respect to nuclear data parameters can be computed efficiently by utilizing the adjoint system of the criticality equation. The propagated nuclear data uncertainty can then be calculated deterministically by the Sandwich rule.

Initially, classical perturbation theory was implemented to the code, which enabled sensitivity analysis of the critical eigenvalue. In this context, covariance matrices from scale 6 were transformed to become compatible with CASMO-4, and the resulting covariance library was connected with the code. Since the cross-section libraries of CASMO-4 do not contain data for the individual capture and scattering reactions, the covariance matrices of the individual subreactions were combined in the covariance library. This work has been reported in detail in [4], and the main conclusions related to the methodology were summarized in this paper. In particular, the sensitivities with respect to total scattering and individual scattering cross-sections cannot be defined in a consistent manner, which leads to some systematic differences in the uncertainty results.

Recently, generalized perturbation theory was added to the code as a new feature, which enables performing sensitivity analysis for responses that can be represented as reaction rate ratios. For each response, the computation of sensitivity profiles with respect to all parameters of interest requires solving one generalized adjoint system. The mathematical background as well as the physical interpretation of the generalized adjoint solutions were reviewed, and practical guidelines were given for modifying a deterministic transport code to solve the generalized adjoint systems needed in sensitivity analysis. The theory for computing the sensitivity profiles was presented both from the perspective of function space analysis and numerical computations.

Numerical results were presented for a lattice physics test problem representing a BWR in hot zero power conditions, and they were compared to the results given by the TSUNAMI-2D sequence in SCALE 6.1. Two-group homogenized cross-sections were considered as responses in the generalized perturbation theory framework. The results were in very good agreement with the thermal responses, whereas in the case of fast responses, the uncertainties given by TSUNAMI-2D were consistently greater. Detailed sensitivity and uncertainty results were presented and analyzed for the homogenized fast and thermal production cross-sections. The differences in the uncertainty results for the fast responses were explained by the incompatible constraints used in computing the scattering uncertainties.

In the future, the work will continue by extending the GPT framework to other responses in addition to two-group homogenized cross-sections with the eventual goal of modifying CASMO-4 to provide uncertainty estimates for all homogenized assembly data, which can then be propagated to coupled neutronics/thermal hydraulics calculations.

Acknowledgment

This work was funded through the Finnish National Research Programme on Nuclear Power Plant Safety 2011–2014, SAFIR2014.