- Email: [email protected]
Generalized equivalence methods for 3D multi-group neutron transport
Annals of Nuclear Energy 112 (2018) 9–16
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Generalized equivalence methods for 3D multi-group neutron transport Guillaume Giudicelli ⇑, Kord Smith, Benoit Forget Massachusetts Institute of Technology, Nuclear Science and Engineering Department, 77, Massachusetts Avenue, Cambridge 02139, United States
a r t i c l e
i n f o
Article history: Received 19 June 2017 Received in revised form 14 September 2017 Accepted 15 September 2017
Keywords: Neutron transport Self-shielding Equivalence Condensation Discontinuity factors Jump conditions
a b s t r a c t The process of generating multi-group cross section data to be used in full core 3D transport models requires not only accurate resonance self-shielding methods, but also some form of equivalence method in order to precisely preserve reaction rates of spectral geometry calculations. This paper presents extensions of the traditional concepts of local reaction rate preservation (common in discontinuity factor, SPH, and BBH methods), to derive a new state-of-the-art transport equivalence method that incorporates angular flux jump conditions that provide polar angle neutron current preservation. This method is tested on numerous fixed-source pin-cell problems by condensing fine energy resonance fluxes and cross sections. The method is demonstrated to precisely preserve all spectral geometry multi-group reaction rates as well as polar angle neutron currents for a wide range of cross section resonance heights, fuel pin diameters, coolant densities, and group energy widths. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction Neutron transport calculations of nuclear reactor cores require accurate methods for treating the complexities in space and energy that arise from the resolved resonances in nuclear cross sections of reactor fuel and structural materials. Treatment of resolved resonance effects is implicit in continuous-energy Monte Carlo methods that represent resonances with hundreds of thousands of energy points. Unfortunately, full-core, steady state reactor calculations with Monte Carlo require tens of billions of neutron histories (Miao et al., 2016) to converge local fluxes and pin power distributions for cores that have near-unity dominance ratios, which is common in commercial light water reactors (LWRs). In addition, Monte Carlo methods are not easily extendable to transient reactor safety analysis where transients must be followed for 100 s or 1000 s of seconds. For these reasons, deterministic neutron transport (or diffusion) methods are likely to be used for performing the bulk of nuclear reactor analysis for years into the future. In deterministic neutron transport methods, the fine-energy detail of Monte Carlo is not practical because of the overwhelming amount of data needed to represent group-to-group transfer cross sections. Instead, deterministic neutron transport models employ some form of approximate resonance treatment that enables direct computation of
⇑ Corresponding author. E-mail address: [email protected] (G. Giudicelli). https://doi.org/10.1016/j.anucene.2017.09.024 0306-4549/Ó 2017 Elsevier Ltd. All rights reserved.
multi-group cross sections for core analysis using only 50 to 500 neutron energy groups. Historically, multi-group cross sections for core analysis have been generated for representative spectral geometries (e.g., 1D pin cells or 2D fuel assemblies). In recent years, even full-core spectral geometries have been used in some core analysis systems. However, the common trend among all deterministic methods is that simplified resonance approximations (e.g., equivalence theory, sub-group, equivalent self-shielding method, etc.) are made in the process of producing multi-group cross sections for core transport models. In order to avoid such simplifying assumptions, some researchers have recently employed direct full-core Monte Carlo calculations and massive cross section tallies to directly compute multi-group cross sections for downstream deterministic transport calculations (Boyd, 2017; Gunow et al., 2017). Such approaches are feasible because Monte Carlo cross section tallies converge much more rapidly than local fluxes and power distributions and because cross section tallies can be agglomerated across many spatial regions through use of sophisticated clustering or machinelearning algorithms (Boyd, 2017). However, none of these methods for generating multi-group cross sections, including Monte Carlo, are capable of directly producing multi-group cross sections that preserve local reaction rates and pin power distributions in 3D core transport models. Even formally exact multi-group cross sections lead to errors in downstream core models (Boyd, 2017; Boyd et al., 2017) unless additional equivalence parameters are introduced to force conservation of local reaction rates.
10
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
2. Equivalence methods The starting point for all equivalence methods is the concept of preservation of reaction rates. If one assumes that a very fine spatial, angular, and energy mesh solution to the neutron transport equation, Eq. (1), is known for the spectral geometry of interest (from either Monte Carlo or from some approximate resonance model), then one can simply define the multi-group transport equation and its corresponding multi-group cross sections as shown in Eqs. (2) and (3).
~; EÞ þ Rt ð~ ~; EÞ ¼ Q ð~ ~r ~ wð~ X r; X r; EÞwð~ r; X r; X; EÞ
ð1Þ
The first term on the left is the streaming of the neutrons, the second the sink rate by absorption and scattering. On the right, Q ð~ r; X; EÞ includes the creation of neutrons by fission, the scattering source and fixed sources. Each term of the neutron transport equation can be integrated over energy groups.
~ ~r ~ w ð~ ~~ ~~ ~ X g r; XÞ þ Rtg ðr; XÞwg ðr; XÞ ¼ Q g ðr; XÞ ~Þ ¼ Rtg ð~ r; X
R E2 E1
~ ; E0 Þ dE0 Rt ð~ r; E0 Þwð~ r; X R E2 0 ~ ; E0 Þ dE wð~ r; X E1
ð2Þ ð3Þ
Each multi-group cross section is then a function of angle. Since hundreds of discrete angles are typically used in reactor core calculations, there is a massive data explosion associated with using multi-group cross sections with a full angular dependency, even if it can be partially mitigated by using polynomial expansions. Consequently, it is common to define cross sections that are angle independent, as shown in Eq. (4), and to use these cross sections in conjunction with Eq. (2) for downstream core analysis.
R E2
Rg ð~ rÞ ¼
E1
dE0 Rð~ r; E0 ÞUð~ r; E0 Þ Rg ¼ R E2 0 Ug dE Uð~ r; E0 Þ E1
ð4Þ
This angular independence of cross section is an approximation, and as a consequence, even a fully converged (in phase space) multi-group calculation will not preserve the desired reaction rates. Equivalence parameters are introduced in order to force this desired preservation of reaction rates. In this process, it is also common to integrate over space in addition to energy and angle over finite spatial volumes (e.g., individual pins, or cladding, or coolant channels), as shown in Eq. (5), and the desired equivalence parameters are now sought to implicitly correct for angularintegration, energy-condensation, and spatial-homogenization.
R R E2
Rg ð~ rÞ ¼
V
dr 0 dE0 Rt ð~ r 0 ; E0 ÞUð~ r 0 ; E0 Þ Rg ¼ R R E2 Ug dr 0 dE0 Uð~ r 0 ; E0 Þ V E1
E1
ð5Þ
The need for practical equivalence parameters has been understood for more than 50 years, and numerous schemes have been proposed and utilized. Bell et al. (1967) proposed a method of introducing artificial anisotropic scattering cross sections (the BHS model) to treat multi-group cross sections with low-order Legendre expansions as the equivalence parameters. While feasible, this method has not been widely adapted in any practical reactor analysis system. For systems that utilize the diffusion approximation for 3D core analysis, Koebke’s heterogeneity factors (Koebke, 1978) and Smith’s extension, discontinuity factors (Smith, 1984) have been widely used equivalence methods for nodal diffusion reactor analysis of commercial LWRs. Discontinuity factors are simply jump conditions on the scalar flux that are applied when the diffusion equations are solved, and the appropriate values of the discontinuity factors are defined to preserve the surface-integrated net neutron current between spatially homogenized and/or energy
collapsed regions. Discontinuity factors are computed directly from lattice spectral geometries used in nodal diffusion core models. In addition, discontinuity factors can be computed for pin-cell spatial regions and used in pin-cell-based 2D or 3D full core reactor models. Discontinuity factors must be saved for each surface of each homogenized/collapsed region, in each of the neutron energy groups. Unfortunately, discontinuity factors are only directly applicable to diffusion theory and not to transport theory equivalence. They have however been derived for simplified P3 equivalence (Kozlowski et al., 2011) and transport acceleration (Grassi, 2007). The most widely used of the transport equivalence methods is Kavenoky’s (Kavenoky, 1978) super-homogenization and its later generalization by Hebert and Benoist (1991), often referred to as the SPH method. The SPH method introduces a factor, f i g, multiplying the multi-group cross sections, and this factor is chosen to conserve the reaction rates in a spatial region i and energy group g. As shown in Eq. (6), the flux computed by the transport solver, USPH ig , is not the same as the group flux from condensation of the reference flux, but rather, it is the flux that is obtained when the transport equation is solved when employing f ig Rig rather than Rig as the group cross section.
ðf ig Rig ÞUSPH ¼ Rig Uig
ð6Þ
In the SPH method, f ig is iteratively determined by solving the multi-group problem over and over again until converged values of f ig are obtained for each spatial region of the spectral geometry. Because the same SPH factor is used on all cross sections within group g, this method preserves up- and down-scatter rates, as well as total reaction rates. In practice, the larger the correction, f ig , the more the SPH flux changes from the reference group flux, and iterative convergence is not guaranteed. The SPH method is very general and it can be used for condensation, for spatial homogenization, or even for employing lower-order transport models. The SPH method is transport-solver-agnostic, and has virtually no memory-cost, since equivalence factors can be absorbed directly into the multi-group cross sections. In an extension to SPH and Discontinuity factors, Sanchez (2009) has defined a class of methods he refers to as black box homogenization (BBH), in which two discontinuity factors are used on each surface of each homogenized/condensation region in order to preserve the neutron partial currents. Sanchez’s implementation of BBH for pin-cell homogenization involves determining the transmission matrix of the system and the escape vector of the internal source, for each pin in the assembly, for each neuron energy group. The BBH class of methods is very general, and it allows preservation of reaction rates in any arbitrary spatial region and energy group of the spectral geometry. In this paper we introduce a new method, which uses an extension of the discontinuity factor concept by applying jump conditions on the angular flux. Jump conditions are explicitly treated as a function of the polar angle neutron fluxes. These equivalence parameters are defined to preserve neutron flow in the polar direction something that neither SPH nor BBH currently provides. The new equivalence methods preserves 3D transport effects, which are anticipated to enable more accurate full-core deterministic transport calculations. These new equivalence methods are, in principle, applicable to any deterministic resonance cross section method or core transport model, but the methods will be demonstrated here using a simple Single Level Breit-Wigner (SLBW) resonance models, in the low temperature approximation, and 2D Step Characteristic pin-cell transport calculations.
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
11
3. A new equivalence method The accuracy of equivalence methods is normally demonstrated by taking a reference solution for some spectral geometry and then performing spatial homogenization and/or energy collapsing to obtain equivalence parameters. One then evaluates the errors that result from the solution of the desired multi-group transport model employing those equivalence parameters for the same spectral geometry that generated the reference solution. Since it is well known that discontinuity factors, SPH factors, and BBH factors can treat either energy condensation, spatial homogenization, or transport resolution reduction, in this paper we focus on the most difficult of those effects by examining the energy condensation consisting of a reduction from fine-energy resolved resonance transport solutions to coarse-energy multi-group transport models. The spectral geometry chosen is the one used by Gibson et al. (2015), in their analysis of resonance self-shielding methods. Gibson modeled an isolated Single-Level Breit-Wigner (SLBW) resonance using the narrow resonance (NR) model for a single fuel pin surrounded by coolant having energy-independent slowing down sources proportional to the potential scattering cross section in each of the two regions. The pin-cell geometry boundary conditions are reflective, which produces a tremendous azimuthal variation of the neutron flux because of the pin-to-pin flux selfshielding. This simple resonance model, shown in Fig. 1a, enables rapid examination of various resonance heights, geometric and material parameters as well as widths of energy groups. SLBW data has been selected to encompass resonances similar to the large near-thermal resonances of U-238. The group widths chosen cover what is typically used in LWR modeling of such resonances. While this is a particularly simple resonance model, it introduces all of the complexity that challenges equivalence methods. The reference solution to this spectral geometry problem is generated using a 2D MOC solver (step characteristics, with spatial flat sources in both fuel and coolant) with a very fine energy grid to map out the large resonance cross section variation in energy. Typical energy dependence of the reference average scalar flux in the fuel and moderator are displayed in Fig. 1b, and one can observe the extremely large dip in fluxes near the resonance energy; which is typical of the classic resonance self-shielding effect. Cross sections are condensed over the width of the resonance group and the one-group analogue of the fine-energy-mesh problem is solved using the energy-collapsed cross section data. Equivalence parameters are sought so that the group reaction rate and flux in both the fuel and the coolant match the reference spectral solution. The equivalence methods proposed here introduces discontinuity factors (e.g., jump conditions) that correct fuel pin incoming and outgoing angular fluxes in order to preserve the group flux. Since this method uses the classic flux-volume-weighted and energy-collapsed group cross sections, all reaction rates will necessarily be preserved if the incoming and outgoing angular fluxes are preserved. This method is essentially an extension of the Generalized Equivalence Theory to transport-theory equivalence. Since this method by definition also conserves all regions’ partial currents, it could also be categorized as falling into the general category of the Black Box Homogenization methods.
3.1. Defining the jump conditions In the method we are using to solve the transport equation, the Method Of Characteristics (MOC), the angular fluxes on which one can apply a jump condition are the track angular fluxes. The most
Fig. 1. Energy dependence in the reference solution. The fuel flux dips where the fuel absorption cross section peaks.
simple idea would be to have a discontinuity factor per track, to make all the track angular fluxes equal to the reference values.
wDF k;a ¼ DF k;a wk;a DF k;a ¼
wref k;a;g
ð7Þ ð8Þ
wk;a
with k the track index, a the angle index and g meaning that it’s averaged over group G, as in
Z
wref k;a;g ¼
Eg
0
0 dE wref k;a ðE Þ=DEg
ð9Þ
ref means the quantity has been computed in the reference solution previously obtained. This would ensure that all track angular fluxes are corrected/ preserved. Preserving them all means getting the same contribution of each track to the scalar flux as the reference solution. This means having the reference group flux as the scalar flux, and so the group reaction rates should be preserved too since we use the group cross-sections. However this would have an even higher memory cost than having an angular group cross section, because there are many tracks per angle and we have one factor per track. So we’ll introduce instead a partial current based jump-condition:
12
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
DF aijg ¼
J ref aijg bJ aijg
ð10Þ
J ref aijg is the contribution of all tracks of angle a to the partial current from region i to j and in group g. This quantity comes from the reference solution. bJ aijg is similar but provided by the simulation using the equivalence method. Fig. 3a illustrates how to obtain these partial currents in the reference calculation when using MOC. Storing these currents has a similar memory cost as storing angular cross sections. By multiplying every track angular flux at the interface between regions by these quantities, we rigorously preserve the partial currents, as the denominator of the discontinuity factor cancels out during the tallying of the partial currents. Preserving the partial currents preserves the reaction rates, a simple neutron balance over the cell guarantees it. These ‘‘angular partial currents” or contribution to the total partial current at the boundary between cell i and j by tracks of angle a, in group g, can be written and approximated in the MOC calculation as follows:
J ref aijg ¼
X
! !
~ k n a wref xa x kaijg
k
ð11Þ
X
xa xk sinðhÞwref kaijg
k
wref kaijg
is the angular flux on track k in group g, at the boundary !
between region i and j. n is the normal to the boundary, heading ! !
!
!
!
outwards. n a ¼ sinðhÞ n X, with X the track direction in the 2D ~ k is the area (in 3D)/length (in 2D) of plane, as shown in Fig. 2. x the surface that is around the track at the intersection with the track. Fig. 2 illustrates these quantities in the case of a track crossing a flat cell boundary. If the track spacing is very small, !
!
1 ~ k cosðh cosðhk Þxk , which is easier to compute. The contrin Xx Þ k
bution of each track to the angular partial current is then similar to the contribution of each track to the scalar flux, written in Eq. 13, which facilitates implementation in a MOC transport solver. These jump conditions are then applied on the track angular fluxes, as the track crosses the surface between region i and j:
^ kaijg wkaijg ¼ DF aijg w
ð12Þ
This process is illustrated in Fig. 3b. To further support the explanation, pseudo-code describing the use of jump conditions in an MOC solver is laid out in Algorithm 1. Steps bolded are added to a regular MOC solver. The reaction rates are defined by multiplying the scalar fluxes by the group cross sections. Each segment of a track in a region contributes to the scalar flux. In the MOC sweep algorithm for every track the scalar flux is incremented by:
Fig. 2. Definition of variables for track k’s contribution to the partial current when crossing a surface, here shown vertical and flat (Boyd et al., 2014).
^ ig þ ¼ U
4p xa xk sinðhÞDw^ kaig Ai
ð13Þ
with Ai the area of the region, xa the angle weight in the angular quadrature, xk the track weight, sinðhÞ the projection on the 2D ^ k;a;i;g ¼ w ^ outgoing w ^ incoming , the variation of the track plane and Dw k;a;i;g
k;a;i;g
angular flux over the segment. To account for the use of discontinuity factors on the track incoming and outgoing angular fluxes when tallying the scalar flux, the variation of flux along the track is redefined as:
^ kaig ¼ DF outgoing w ^ outgoing DF incoming w ^ incoming Dw aig aig kaig kaig
ð14Þ
incoming means at the location where the track enters region i, and outgoing means at the location where the track exits region i. Algorithm 1 Pseudo-code for MOC algorithm with jump conditions. Initialize scalar fluxes while scalar fluxes not converged do for loop over azimuthal angles do for loop over tracks do for loop over segments do for loop over energy groups do for loop over polar angles do Compute angular flux change along segment Tally contribution to angular partial current Compute discontinuity factor from Eq. (10) Apply discontinuity factor on outgoing flux Update angular flux change along segment Increment scalar flux Update track outgoing flux end for end for end for end for end for Set incoming fluxes with boundary conditions Update keff and FSR sources Update previous iteration angular partial currents endwhile
In the definition of the discontinuity factors that we proposed, the numerator is given by the reference calculation, and the denominator is accumulated at every sweep. However, in the transport sweep, we don’t have access to the entire sum that makes up the denominator as we are doing the sweep, which is only known after all tracks have been swept. Instead, we propose to use the angular partial current (the DF denominator) from the previous iteration. During the first iteration, the boundary angular fluxes are not converged, so the denominator is very far from converged, and the discontinuity factor is very far from one, which leads to a few possible observations: If one uses discontinuity factors early in the calculation, it could help speed up convergence since part of the reference solution is inserted. Unless we start with converged boundary fluxes, or a good guess, using condensed boundary fluxes for example, the discontinuity factor cannot be determined in a single transport sweep.
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
13
Fig. 3. Graphical representation of how to obtain the jump conditions in MOC and how to apply them.
4. Test case We implemented the proposed new method in a 2D MOC toy code. This 2D code was compared with OpenMOC (Boyd et al., 2014) to verify accuracy. The test-case is a square 2D pin-cell with a pitch of 1.26 cm and a pin radius of 0.4 cm unless stated otherwise. There is resonant absorption in the fuel, but resonance scattering is ignored in the initial development. The energy group width is 4 eV unless mentioned otherwise, and the resonance peak is located at 1000 eV in the center of the energy group. The cross section data is summarized in Table 1. The bare fuel pin model needed 128 azimuthal angles (Legendre quadrature), a track spacing of 0.001 cm and 3 polar angles (TY quadrature (Yamamoto et al., 2007)) to get results converged to a few pcms. For a lattice, the number of tracks and azimuthal angles needed is less as the tracks sample the pins multiple times when they cross the geometry. 5. Implementation and results At this point, the proposed method uses a discontinuity factor (DF) per angle and per region. This strategy has the same memory requirement as using angular cross sections, which would be preferred to using DFs. The numerator and denominator of these discontinuity factors can be averaged over groups of angles and then the same factor can be used for all angles in the group. This lowers significantly the memory cost. In Table 2, we can see the error on the fuel absorption reaction rates without equivalence (column 2) and with equivalence (column 3 to 9). On the left part of the table, we have applied discontinuity factors on the outgoing track angular fluxes only, and on the right side on the incoming and outgoing track angular fluxes. Then, in each of these two parts, we have allowed the discontinuity factors to depend either on: -azimuthal and polar angles and regions or -polar angles and regions only or -regions only (no angular dependence). When not depending on an angle set, we simply average the discontinuity factors over the angle set using the rele-
Table 1 Densities and non-resonant cross section data for the test case. Material
q(g/cc)
Nuclide
rsc: (b)
rabs: ðbÞ
UO2 H2O
10.4 0.7
U-238 O-16 H-1
11.29 3.888 20.47
0 0 0
vant quadrature. We can see that applying discontinuity factors on the incoming and outgoing fluxes yields the best results. When applying jump conditions depending on polar angles only, we do not lose accuracy compared with having azimuthally and polar angle dependent discontinuity factors. The discontinuity factors can have a polar and an azimuthal dependency, but a factor per surface independent of angle is enough to preserve partial currents. The polar dependence would be a unique feature compared to traditional equivalence methods because it accounts for axial effects. The reader should note that the errors mentioned are fractional errors on the fuel absorption rate for a single resonance in a single fuel pin. When scaling up to a full reactor problem, the maximum error coming from the flux separability approximation on U 238 absorption rates is closer to 2% and the mean error on eigenvalue closer to 400 pcm (Boyd, 2017). After the number of angles, the second factor of great memory costs is the number of regions. To decrease the memory requirements, we could make discontinuity factors a function of their neighborhood (Herrero et al., 1997), possibly with an artificial nodal network (Kozlowski et al., 2004), we could use clustering algorithms for discontinuity factors, like for the group crosssections (Boyd, 2017), or we could seek ways to propagate the factors inside a fuel pin from data gathered over a single pin region.
5.1. Parametric study At this point we want to show that this new method can be applied to a variety of cases derived from the resonant pin-cell problem. For that we perform a parametric study over pin radius, moderator density, presence of resonant scattering in the fuel and finally energy group width. Given the previous section’s results, we will only study polar angle dependent jump conditions applied on the incoming and outgoing track angular fluxes. So far, they are the best compromise between memory usage and accuracy of the reaction rates, with some potential for 3D equivalence due to the polar angle dependence.
5.1.1. Varying the fuel pin radius With increasing pin radius, the amount of self-shielding increases, so errors on the fuel absorption cross section should be larger. Table 3 shows results for varying pin radii. Without an equivalence method, the error first gets larger as expected without the correction, then reduces. This reduction is due to the compet-
14
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
Table 2 Fuel absorption reaction rate fractional errors (pcm) for different dependencies of the partial currents-based DFs. The error on the moderator flux is of the same order of magnitude. Discontinuity factors need to be applied on incoming and outgoing angular fluxes to preserve reaction rates. An angular dependence of DFs is not required to achieve conservation of reaction rates. DF used on:
No DF
Outgoing angular fluxes
DF dependencies:
No DF
Surface polar azimuthal
Surface polar
Surface
Incoming and outgoing angular fluxes Surface polar azimuthal
Surface polar
Surface
25 238 828 2206
6 23 47 29
3 15 32 38
9 47 115 153
1 1 1 0
0 0 0 0
0 0 0 0
rpeak ðbÞ 1e2 1e3 1e4 1e5
ing effect of the slowing down source in the fuel becoming more and more important as the fuel radius increases. For the discontinuity factors, first let’s observe that they are inducing much more correction on the pin-outgoing track angular fluxes than on the incoming and that the factors are usually above one. As previously stated, this is because there is insufficient selfshielding in the fuel when using the group cross-sections without their spatial and angular dependence. In order to reduce the group flux, the outgoing current has to be increased, which is why the factors are mostly above one. The discontinuity factors follow mostly the same trends as the error without equivalence. Then, we can see that the steepest and flattest polar angles exhibit opposite trends. The most horizontal polar angle discontinuity factors increases with pin radius, whereas the most vertical decreases with pin radius. This is driven by the source term in the fuel flux. The first/steepest polar angle is the one with the longest path in the fuel, so the contribution of the source is the largest for this angle. On the other hand, the length in the fuel along a track with the third polar angle is the shortest, so the contribution of the source compared to the contribution of the incoming flux is not as big. As mentioned above, the bigger the radius, the more the source in the fuel contributes to the fuel absorption reaction rate, so the lower the required correction. The correction introduced by the jump conditions is much larger for smaller pin radii. For larger pin radii, the distance ‘‘traveled” by the tracks in the moderator between exiting and entering gets smaller, meaning the effect of applying the outgoing jump condition is less dampened by the coolant before re-entering the fuel pin. These two effects could each introduce numerical instability. By numerical instability, we mean that the track angular fluxes could diverge because of consecutive multiplications by the discontinuity factors. As shown in this table, the method has reached convergence in all cases and reduced the error on reaction rates and group fluxes by orders of magnitude. 5.1.2. Varying the moderator density In order to show the applicability of our method for both boiling and pressurized water reactors, at zero power as well as full power, we vary the moderator density and the fuel pin radius simultaneously. Results shown in Table 4 indicate that the accuracy of our
method is not affected. As expected, the higher the water density, the stronger the source in the water and the higher the error in the absorption rate due to the error in self-shielding when using the group-cross section without equivalence factors. 5.1.3. Introducing resonant scattering in the fuel So far we have only modeled resonant absorption in the fuel. A more realistic model would include resonant scattering. This could create issues with the equivalence, as now the error comes both from the condensation of the absorption and scattering crosssections. In order to capture the up-scatter and down-scatter rates of neutrons properly it is important to get the scattering rates right. Our method preserves group fluxes, so it preserves all reaction rates indiscriminately. To further support this claim, we present in Table 5 results with and without resonant scattering in the fuel. In order to introduce scattering without changing the absorption cross section energy profile, we choose to adapt the resonance widths of scattering and absorption. The total width is kept at 0.025 eV, while the absorption width is 0.015 eV and the scattering width 0.01 eV. The resonant scattering energy dependence is then exactly the same as in Fig. 1a. We first observe that the introduction of scattering diminishes the absorption group cross section. In the narrow resonance model, when neutrons scatter, they scatter out of the resonance and can’t be absorbed anymore. Even though there is less absorption going on, the error when not using equivalence methods increases in the presence of scattering. Resonant scattering in the narrow resonance model acts very similarly as absorption, as only potential scattering contributes to the source term. So the self-shielding effectively increases with the presence of scattering, and the error on absorption rates gets larger. Finally, as the self shielding increases, the equivalence factors are getting larger to correct the increasing error on the group fluxes and reaction rates. 5.1.4. Varying energy group widths Because the effect of resonant absorption on the absorption group cross section strongly depends on the energy group width, we need to make sure that our method is applicable for the wide range of energy group widths that is used in current deterministic transport solvers. In Table 6, the method is tested for groups of
Table 3 Fractional errors on fuel absorption reaction rate and polar discontinuity factors for various radii of a fuel pin with a 104 b peak SLBW absorption-only resonance in the fuel. DFs on the incoming angular fluxes only perform minor corrections, and DFs on the outgoing angular fluxes are larger for cases with larger errors without equivalence. Polar angle 1 is the steepest in the quadrature. Pin radius (cm)
Error on fuel absorption rate (pcm)
Polar incoming DF
Polar outgoing DF
No DF
DF
1
2
3
1
2
3
0.2 0.3 0.4 0.5 0.6
859 910 829 621 261
<0.1 <0.1 <0.1 <0.1 <0.1
1.00001 1.00002 1.00009 1.00027 1.00067
1.00011 1.00021 1.00032 1.00031 1.00051
1.00024 1.00019 1.00010 0.99995 1.00017
1.14989 1.13096 1.10441 1.07431 1.04102
1.02155 1.02266 1.02103 1.01684 1.01028
0.99671 0.99827 0.99999 1.00132 1.00167
15
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
Table 4 Fractional errors on fuel absorption reaction rate for various radii of a fuel pin at various water densities with a 104 b peak SLBW absorption-only resonance in the fuel. The pin-topin pitch is still 1.26 cm. DFs allow for exact preservation of reaction rates in all cases, representative of LWR operating conditions. Fractional error on fuel absorption rate (pcm) 0.2 g/cc Pin rad.
0.7 g/cc
1 g/cc
No DF
DF
No DF
DF
No DF
DF
643 580 401 209 50
<0.1 <0.1 <0.1 <0.1 <0.1
859 910 829 621 261
<0.1 <0.1 <0.1 <0.1 <0.1
896 979 936 754 363
<0.1 <0.1 <0.1 <0.1 <0.1
qmod 0.2 0.3 0.4 0.5 0.6
Table 5 Fractional errors on fuel absorption reaction rate (in pcm), absorption group cross sections and polar outgoing discontinuity factors with resonant absorption and without/with resonant scattering. The introduction of scattering shields the absorption but also increases the self-shielding error. Discontinuity factors increase to provide additional correction and allow for exact preservation of reaction rates. Without resonant scattering
With resonant scattering
rpeak
rabs;g
Fuel abs. error
Polar outgoing DF
(b)
(b)
No DF
DF
Angle 1
2
1e1 1e2 1e3 1e4 1e5
.09 .63 2.4 8.3 30
0.4 25 238 828 2206
<0.1 <0.1 <0.1 <0.1 <0.1
1.00029 1.00703 1.03362 1.10441 1.23649
1.00002 1.00084 1.00567 1.02103 1.08287
rabs;g
Fuel abs. error
Polar outgoing DF
3
(b)
No DF
DF
Angle 1
2
3
0.99999 0.99983 0.99954 0.99999 1.01376
.09 .55 2.0 6.8 25
3.2 60 357 1168 3059
<0.1 <0.1 <0.1 <0.1 0.9
1.00205 1.01590 1.06050 1.19302 1.49169
1.00014 1.00188 1.00881 1.03188 1.13024
0.99995 0.99964 0.99913 0.99917 1.01872
Table 6 Fractional errors on fuel absorption reaction rate and polar discontinuity factors for various energy group-widths (eV) with a 104 barns peak SLBW absorption-only resonance in the fuel. The smaller the energy group width, the higher the error on the absorption reaction rate, and the larger the discontinuity factors. Group width (eV)
20 10 4 2 1 0.5
Fuel abs. error (pcm)
Polar incoming DF
No DF
DF
1
2
3
1
2
3
174 345 829 1535 2504 2737
<0.1 <0.1 <0.1 <0.1 <0.1 <0.1
1.00002 1.00005 1.00009 1.00010 1.00005 0.99999
1.00008 1.00015 1.00032 1.00052 1.00070 1.00055
1.00003 1.00005 1.00010 1.00013 1.00017 1.00048
1.02241 1.04416 1.10441 1.18303 1.24798 1.17955
1.00381 1.00781 1.02103 1.04687 1.10906 1.21131
0.99972 0.99957 0.99999 1.00373 1.02396 1.11074
widths ranging from 20 to 0.5 eV, centered on a single SLBW resonance. We can see that the condensation error without an equivalence method gets larger as the group width decreases. Simultaneously, the discontinuity factors increase mostly monotonously to correct the partial currents.
6. Conclusion and future work Cross sections for multi-group transport calculations are traditionally generated using the scalar flux, thus neglecting the angular dependence of group cross sections. This approximation leads to a large systematic bias on pin power distributions in light water reactor systems, and an equivalence method has to be used to preserve reaction rates. This paper introduces a new state-of-the-art equivalence method that incorporates an extended level of transport solution conservation; namely direct polar angle current preservation. It is based on jump conditions applied on the incoming and outgoing track angular fluxes at the fuel pin surface. These jump conditions are chosen to be polar angle dependent to enforce polar angle partial currents preservation. This feature is expected to be important for 3D MOC transport calculations. The method is tested on a 2D fixed source pin-cell problem by condensing fine energy resonance fluxes and resonant cross sec-
Polar outgoing DF
tions in a parametric study on pin radius, moderator density, presence of resonant scattering and energy group width. The parameter ranges were chosen to encompass LWR operating conditions, and common group widths in deterministic multi-group transport. The method demonstrates precise preservation of all reaction rates in all cases. This equivalence method has already been extended to the case of spatial representations with multiple fuel rings commonly used in fuel depletion and fuel temperature distribution modeling for LWRs and also showed exact preservation. While this article assumed reference partial currents were known in order to derive equivalence parameters, future work will focus on alternative ways to obtain the jump conditions. This work will also soon be extended to the case of Monte Carlo-generated group cross sections so that one can obtain accurate equivalence parameters directly from partially converged, full-core, Monte Carlo calculations; thereby avoiding all approximations inherent in traditional multi-group cross section generation methods.
Acknowledgements The research work presented in this paper was supported by the US Department of Energy Nuclear Energy University Program contract DE-NE0008578.
16
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
References Bell, G.I., Hansen, G.E., Sandmeier, H.A., 1967. Multitable treatments of anisotropic scattering in SN multigroup transport calculations. Nucl. Sci. Eng. 379, 28. Boyd, W., 2017. Reactor Agnostic Multi-Group Cross Section Generation for FineMesh Deterministic Neutron Transport Simulations (Ph.D. thesis). Massachusetts Institute of Technology. Boyd, W., Gibson, N., Forget, B., Smith, K., 2017. An analysis of condensation errors in multi-group cross section generation for fine mesh neutron transport calculations. Submitted to Annals of Nuclear Energy. Boyd, W., Shaner, S., Li, L., Forget, B., Smith, K., 2014. The openmoc method of characteristics neutral particle transport code. Ann. Nucl. Energy 68, 43. Gibson, N., Smith, K., Forget, B., 2015. Simple benchmark for evaluating selfshielding models. In: Proceedings of M&C 2015, Nashville, Tennessee. Grassi, G., 2007. A nonlinear space-angle multigrid acceleration for the method of characteristics in unstructured meshes. Nucl. Sci. Eng. 155 (2), 208–222. Gunow, G., Shaner, S., Boyd, W., Forget, B., Smith, K., 2017. Accuracy and performance of 3D MOC for full-core PWR problems. In: Proceedings of M&C2017, Jeju Korea. Hebert, A., Benoist, P., 1991. A consistent technique for the global homogenization of a pressurized water reactor assembly. Nucl. Sci. Eng. 109 (4), 360–372. Herrero, J., Garcia-Herranz, N., Ahnert, C., 1997. Neighborhood corrected interface discontinuity factors for multi-group pin-by-pin diffusion calculations for LWR. Ann. Nucl. Energy 24 (6), 477–488.
Kavenoky, A., 1978. The SPH homogenization method. In: IAEA Proceedings of Specialists Meeting on Homogenization Method in Reactor Physics, Lugano, Switzerland. Koebke, K., 1978. A new approach to homogenization and group condensation. In: IAEA Proceedings of Specialists Meeting on Homogenization Method in Reactor Physics, Lugano, Switzerland. Kozlowski, T., Lee, D., Downar, T.J., 2004. The use of artificial neural network for online prediction of pin-cell discontinuity factors in PARCS. In: PHYSOR 2004 Proceedings. Kozlowski, T., Xu, Y., Downar, T.J., Lee, D., 2011. Cell homogenization method for pin-by-pin neutron transport calculations. Nucl. Sci. Eng. 169, 1–18. Miao, J., Forget, B., Smith, K., 2016. Analysis of correlations and their impact on convergence rates in Monte Carlo eigenvalue simulations. Ann. Nucl. Energy 92, 81–95. Sanchez, R., 2009. Assembly homogenization techniques for core calculations. Prog. Nucl. Energy 51, 14–31. Smith, K., 1984. Nodal method storage reduction by non-linear iteration. Trans. ANS 44, 265. Yamamoto, A., Tabuchi, M., Sugimura, N., Ushio, T., Mori, M., 2007. Derivation of optimal polar angle quadrature set for the method of characteristics based on approximation error of the Bickley functions. J. Nucl. Sci. Technol. 44 (2), 129– 136.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Generalized equivalence methods for 3D multi-group neutron transport Guillaume Giudicelli ⇑, Kord Smith, Benoit Forget Massachusetts Institute of Technology, Nuclear Science and Engineering Department, 77, Massachusetts Avenue, Cambridge 02139, United States
a r t i c l e
i n f o
Article history: Received 19 June 2017 Received in revised form 14 September 2017 Accepted 15 September 2017
Keywords: Neutron transport Self-shielding Equivalence Condensation Discontinuity factors Jump conditions
a b s t r a c t The process of generating multi-group cross section data to be used in full core 3D transport models requires not only accurate resonance self-shielding methods, but also some form of equivalence method in order to precisely preserve reaction rates of spectral geometry calculations. This paper presents extensions of the traditional concepts of local reaction rate preservation (common in discontinuity factor, SPH, and BBH methods), to derive a new state-of-the-art transport equivalence method that incorporates angular flux jump conditions that provide polar angle neutron current preservation. This method is tested on numerous fixed-source pin-cell problems by condensing fine energy resonance fluxes and cross sections. The method is demonstrated to precisely preserve all spectral geometry multi-group reaction rates as well as polar angle neutron currents for a wide range of cross section resonance heights, fuel pin diameters, coolant densities, and group energy widths. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction Neutron transport calculations of nuclear reactor cores require accurate methods for treating the complexities in space and energy that arise from the resolved resonances in nuclear cross sections of reactor fuel and structural materials. Treatment of resolved resonance effects is implicit in continuous-energy Monte Carlo methods that represent resonances with hundreds of thousands of energy points. Unfortunately, full-core, steady state reactor calculations with Monte Carlo require tens of billions of neutron histories (Miao et al., 2016) to converge local fluxes and pin power distributions for cores that have near-unity dominance ratios, which is common in commercial light water reactors (LWRs). In addition, Monte Carlo methods are not easily extendable to transient reactor safety analysis where transients must be followed for 100 s or 1000 s of seconds. For these reasons, deterministic neutron transport (or diffusion) methods are likely to be used for performing the bulk of nuclear reactor analysis for years into the future. In deterministic neutron transport methods, the fine-energy detail of Monte Carlo is not practical because of the overwhelming amount of data needed to represent group-to-group transfer cross sections. Instead, deterministic neutron transport models employ some form of approximate resonance treatment that enables direct computation of
⇑ Corresponding author. E-mail address: [email protected] (G. Giudicelli). https://doi.org/10.1016/j.anucene.2017.09.024 0306-4549/Ó 2017 Elsevier Ltd. All rights reserved.
multi-group cross sections for core analysis using only 50 to 500 neutron energy groups. Historically, multi-group cross sections for core analysis have been generated for representative spectral geometries (e.g., 1D pin cells or 2D fuel assemblies). In recent years, even full-core spectral geometries have been used in some core analysis systems. However, the common trend among all deterministic methods is that simplified resonance approximations (e.g., equivalence theory, sub-group, equivalent self-shielding method, etc.) are made in the process of producing multi-group cross sections for core transport models. In order to avoid such simplifying assumptions, some researchers have recently employed direct full-core Monte Carlo calculations and massive cross section tallies to directly compute multi-group cross sections for downstream deterministic transport calculations (Boyd, 2017; Gunow et al., 2017). Such approaches are feasible because Monte Carlo cross section tallies converge much more rapidly than local fluxes and power distributions and because cross section tallies can be agglomerated across many spatial regions through use of sophisticated clustering or machinelearning algorithms (Boyd, 2017). However, none of these methods for generating multi-group cross sections, including Monte Carlo, are capable of directly producing multi-group cross sections that preserve local reaction rates and pin power distributions in 3D core transport models. Even formally exact multi-group cross sections lead to errors in downstream core models (Boyd, 2017; Boyd et al., 2017) unless additional equivalence parameters are introduced to force conservation of local reaction rates.
10
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
2. Equivalence methods The starting point for all equivalence methods is the concept of preservation of reaction rates. If one assumes that a very fine spatial, angular, and energy mesh solution to the neutron transport equation, Eq. (1), is known for the spectral geometry of interest (from either Monte Carlo or from some approximate resonance model), then one can simply define the multi-group transport equation and its corresponding multi-group cross sections as shown in Eqs. (2) and (3).
~; EÞ þ Rt ð~ ~; EÞ ¼ Q ð~ ~r ~ wð~ X r; X r; EÞwð~ r; X r; X; EÞ
ð1Þ
The first term on the left is the streaming of the neutrons, the second the sink rate by absorption and scattering. On the right, Q ð~ r; X; EÞ includes the creation of neutrons by fission, the scattering source and fixed sources. Each term of the neutron transport equation can be integrated over energy groups.
~ ~r ~ w ð~ ~~ ~~ ~ X g r; XÞ þ Rtg ðr; XÞwg ðr; XÞ ¼ Q g ðr; XÞ ~Þ ¼ Rtg ð~ r; X
R E2 E1
~ ; E0 Þ dE0 Rt ð~ r; E0 Þwð~ r; X R E2 0 ~ ; E0 Þ dE wð~ r; X E1
ð2Þ ð3Þ
Each multi-group cross section is then a function of angle. Since hundreds of discrete angles are typically used in reactor core calculations, there is a massive data explosion associated with using multi-group cross sections with a full angular dependency, even if it can be partially mitigated by using polynomial expansions. Consequently, it is common to define cross sections that are angle independent, as shown in Eq. (4), and to use these cross sections in conjunction with Eq. (2) for downstream core analysis.
R E2
Rg ð~ rÞ ¼
E1
dE0 Rð~ r; E0 ÞUð~ r; E0 Þ Rg ¼ R E2 0 Ug dE Uð~ r; E0 Þ E1
ð4Þ
This angular independence of cross section is an approximation, and as a consequence, even a fully converged (in phase space) multi-group calculation will not preserve the desired reaction rates. Equivalence parameters are introduced in order to force this desired preservation of reaction rates. In this process, it is also common to integrate over space in addition to energy and angle over finite spatial volumes (e.g., individual pins, or cladding, or coolant channels), as shown in Eq. (5), and the desired equivalence parameters are now sought to implicitly correct for angularintegration, energy-condensation, and spatial-homogenization.
R R E2
Rg ð~ rÞ ¼
V
dr 0 dE0 Rt ð~ r 0 ; E0 ÞUð~ r 0 ; E0 Þ Rg ¼ R R E2 Ug dr 0 dE0 Uð~ r 0 ; E0 Þ V E1
E1
ð5Þ
The need for practical equivalence parameters has been understood for more than 50 years, and numerous schemes have been proposed and utilized. Bell et al. (1967) proposed a method of introducing artificial anisotropic scattering cross sections (the BHS model) to treat multi-group cross sections with low-order Legendre expansions as the equivalence parameters. While feasible, this method has not been widely adapted in any practical reactor analysis system. For systems that utilize the diffusion approximation for 3D core analysis, Koebke’s heterogeneity factors (Koebke, 1978) and Smith’s extension, discontinuity factors (Smith, 1984) have been widely used equivalence methods for nodal diffusion reactor analysis of commercial LWRs. Discontinuity factors are simply jump conditions on the scalar flux that are applied when the diffusion equations are solved, and the appropriate values of the discontinuity factors are defined to preserve the surface-integrated net neutron current between spatially homogenized and/or energy
collapsed regions. Discontinuity factors are computed directly from lattice spectral geometries used in nodal diffusion core models. In addition, discontinuity factors can be computed for pin-cell spatial regions and used in pin-cell-based 2D or 3D full core reactor models. Discontinuity factors must be saved for each surface of each homogenized/collapsed region, in each of the neutron energy groups. Unfortunately, discontinuity factors are only directly applicable to diffusion theory and not to transport theory equivalence. They have however been derived for simplified P3 equivalence (Kozlowski et al., 2011) and transport acceleration (Grassi, 2007). The most widely used of the transport equivalence methods is Kavenoky’s (Kavenoky, 1978) super-homogenization and its later generalization by Hebert and Benoist (1991), often referred to as the SPH method. The SPH method introduces a factor, f i g, multiplying the multi-group cross sections, and this factor is chosen to conserve the reaction rates in a spatial region i and energy group g. As shown in Eq. (6), the flux computed by the transport solver, USPH ig , is not the same as the group flux from condensation of the reference flux, but rather, it is the flux that is obtained when the transport equation is solved when employing f ig Rig rather than Rig as the group cross section.
ðf ig Rig ÞUSPH ¼ Rig Uig
ð6Þ
In the SPH method, f ig is iteratively determined by solving the multi-group problem over and over again until converged values of f ig are obtained for each spatial region of the spectral geometry. Because the same SPH factor is used on all cross sections within group g, this method preserves up- and down-scatter rates, as well as total reaction rates. In practice, the larger the correction, f ig , the more the SPH flux changes from the reference group flux, and iterative convergence is not guaranteed. The SPH method is very general and it can be used for condensation, for spatial homogenization, or even for employing lower-order transport models. The SPH method is transport-solver-agnostic, and has virtually no memory-cost, since equivalence factors can be absorbed directly into the multi-group cross sections. In an extension to SPH and Discontinuity factors, Sanchez (2009) has defined a class of methods he refers to as black box homogenization (BBH), in which two discontinuity factors are used on each surface of each homogenized/condensation region in order to preserve the neutron partial currents. Sanchez’s implementation of BBH for pin-cell homogenization involves determining the transmission matrix of the system and the escape vector of the internal source, for each pin in the assembly, for each neuron energy group. The BBH class of methods is very general, and it allows preservation of reaction rates in any arbitrary spatial region and energy group of the spectral geometry. In this paper we introduce a new method, which uses an extension of the discontinuity factor concept by applying jump conditions on the angular flux. Jump conditions are explicitly treated as a function of the polar angle neutron fluxes. These equivalence parameters are defined to preserve neutron flow in the polar direction something that neither SPH nor BBH currently provides. The new equivalence methods preserves 3D transport effects, which are anticipated to enable more accurate full-core deterministic transport calculations. These new equivalence methods are, in principle, applicable to any deterministic resonance cross section method or core transport model, but the methods will be demonstrated here using a simple Single Level Breit-Wigner (SLBW) resonance models, in the low temperature approximation, and 2D Step Characteristic pin-cell transport calculations.
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
11
3. A new equivalence method The accuracy of equivalence methods is normally demonstrated by taking a reference solution for some spectral geometry and then performing spatial homogenization and/or energy collapsing to obtain equivalence parameters. One then evaluates the errors that result from the solution of the desired multi-group transport model employing those equivalence parameters for the same spectral geometry that generated the reference solution. Since it is well known that discontinuity factors, SPH factors, and BBH factors can treat either energy condensation, spatial homogenization, or transport resolution reduction, in this paper we focus on the most difficult of those effects by examining the energy condensation consisting of a reduction from fine-energy resolved resonance transport solutions to coarse-energy multi-group transport models. The spectral geometry chosen is the one used by Gibson et al. (2015), in their analysis of resonance self-shielding methods. Gibson modeled an isolated Single-Level Breit-Wigner (SLBW) resonance using the narrow resonance (NR) model for a single fuel pin surrounded by coolant having energy-independent slowing down sources proportional to the potential scattering cross section in each of the two regions. The pin-cell geometry boundary conditions are reflective, which produces a tremendous azimuthal variation of the neutron flux because of the pin-to-pin flux selfshielding. This simple resonance model, shown in Fig. 1a, enables rapid examination of various resonance heights, geometric and material parameters as well as widths of energy groups. SLBW data has been selected to encompass resonances similar to the large near-thermal resonances of U-238. The group widths chosen cover what is typically used in LWR modeling of such resonances. While this is a particularly simple resonance model, it introduces all of the complexity that challenges equivalence methods. The reference solution to this spectral geometry problem is generated using a 2D MOC solver (step characteristics, with spatial flat sources in both fuel and coolant) with a very fine energy grid to map out the large resonance cross section variation in energy. Typical energy dependence of the reference average scalar flux in the fuel and moderator are displayed in Fig. 1b, and one can observe the extremely large dip in fluxes near the resonance energy; which is typical of the classic resonance self-shielding effect. Cross sections are condensed over the width of the resonance group and the one-group analogue of the fine-energy-mesh problem is solved using the energy-collapsed cross section data. Equivalence parameters are sought so that the group reaction rate and flux in both the fuel and the coolant match the reference spectral solution. The equivalence methods proposed here introduces discontinuity factors (e.g., jump conditions) that correct fuel pin incoming and outgoing angular fluxes in order to preserve the group flux. Since this method uses the classic flux-volume-weighted and energy-collapsed group cross sections, all reaction rates will necessarily be preserved if the incoming and outgoing angular fluxes are preserved. This method is essentially an extension of the Generalized Equivalence Theory to transport-theory equivalence. Since this method by definition also conserves all regions’ partial currents, it could also be categorized as falling into the general category of the Black Box Homogenization methods.
3.1. Defining the jump conditions In the method we are using to solve the transport equation, the Method Of Characteristics (MOC), the angular fluxes on which one can apply a jump condition are the track angular fluxes. The most
Fig. 1. Energy dependence in the reference solution. The fuel flux dips where the fuel absorption cross section peaks.
simple idea would be to have a discontinuity factor per track, to make all the track angular fluxes equal to the reference values.
wDF k;a ¼ DF k;a wk;a DF k;a ¼
wref k;a;g
ð7Þ ð8Þ
wk;a
with k the track index, a the angle index and g meaning that it’s averaged over group G, as in
Z
wref k;a;g ¼
Eg
0
0 dE wref k;a ðE Þ=DEg
ð9Þ
ref means the quantity has been computed in the reference solution previously obtained. This would ensure that all track angular fluxes are corrected/ preserved. Preserving them all means getting the same contribution of each track to the scalar flux as the reference solution. This means having the reference group flux as the scalar flux, and so the group reaction rates should be preserved too since we use the group cross-sections. However this would have an even higher memory cost than having an angular group cross section, because there are many tracks per angle and we have one factor per track. So we’ll introduce instead a partial current based jump-condition:
12
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
DF aijg ¼
J ref aijg bJ aijg
ð10Þ
J ref aijg is the contribution of all tracks of angle a to the partial current from region i to j and in group g. This quantity comes from the reference solution. bJ aijg is similar but provided by the simulation using the equivalence method. Fig. 3a illustrates how to obtain these partial currents in the reference calculation when using MOC. Storing these currents has a similar memory cost as storing angular cross sections. By multiplying every track angular flux at the interface between regions by these quantities, we rigorously preserve the partial currents, as the denominator of the discontinuity factor cancels out during the tallying of the partial currents. Preserving the partial currents preserves the reaction rates, a simple neutron balance over the cell guarantees it. These ‘‘angular partial currents” or contribution to the total partial current at the boundary between cell i and j by tracks of angle a, in group g, can be written and approximated in the MOC calculation as follows:
J ref aijg ¼
X
! !
~ k n a wref xa x kaijg
k
ð11Þ
X
xa xk sinðhÞwref kaijg
k
wref kaijg
is the angular flux on track k in group g, at the boundary !
between region i and j. n is the normal to the boundary, heading ! !
!
!
!
outwards. n a ¼ sinðhÞ n X, with X the track direction in the 2D ~ k is the area (in 3D)/length (in 2D) of plane, as shown in Fig. 2. x the surface that is around the track at the intersection with the track. Fig. 2 illustrates these quantities in the case of a track crossing a flat cell boundary. If the track spacing is very small, !
!
1 ~ k cosðh cosðhk Þxk , which is easier to compute. The contrin Xx Þ k
bution of each track to the angular partial current is then similar to the contribution of each track to the scalar flux, written in Eq. 13, which facilitates implementation in a MOC transport solver. These jump conditions are then applied on the track angular fluxes, as the track crosses the surface between region i and j:
^ kaijg wkaijg ¼ DF aijg w
ð12Þ
This process is illustrated in Fig. 3b. To further support the explanation, pseudo-code describing the use of jump conditions in an MOC solver is laid out in Algorithm 1. Steps bolded are added to a regular MOC solver. The reaction rates are defined by multiplying the scalar fluxes by the group cross sections. Each segment of a track in a region contributes to the scalar flux. In the MOC sweep algorithm for every track the scalar flux is incremented by:
Fig. 2. Definition of variables for track k’s contribution to the partial current when crossing a surface, here shown vertical and flat (Boyd et al., 2014).
^ ig þ ¼ U
4p xa xk sinðhÞDw^ kaig Ai
ð13Þ
with Ai the area of the region, xa the angle weight in the angular quadrature, xk the track weight, sinðhÞ the projection on the 2D ^ k;a;i;g ¼ w ^ outgoing w ^ incoming , the variation of the track plane and Dw k;a;i;g
k;a;i;g
angular flux over the segment. To account for the use of discontinuity factors on the track incoming and outgoing angular fluxes when tallying the scalar flux, the variation of flux along the track is redefined as:
^ kaig ¼ DF outgoing w ^ outgoing DF incoming w ^ incoming Dw aig aig kaig kaig
ð14Þ
incoming means at the location where the track enters region i, and outgoing means at the location where the track exits region i. Algorithm 1 Pseudo-code for MOC algorithm with jump conditions. Initialize scalar fluxes while scalar fluxes not converged do for loop over azimuthal angles do for loop over tracks do for loop over segments do for loop over energy groups do for loop over polar angles do Compute angular flux change along segment Tally contribution to angular partial current Compute discontinuity factor from Eq. (10) Apply discontinuity factor on outgoing flux Update angular flux change along segment Increment scalar flux Update track outgoing flux end for end for end for end for end for Set incoming fluxes with boundary conditions Update keff and FSR sources Update previous iteration angular partial currents endwhile
In the definition of the discontinuity factors that we proposed, the numerator is given by the reference calculation, and the denominator is accumulated at every sweep. However, in the transport sweep, we don’t have access to the entire sum that makes up the denominator as we are doing the sweep, which is only known after all tracks have been swept. Instead, we propose to use the angular partial current (the DF denominator) from the previous iteration. During the first iteration, the boundary angular fluxes are not converged, so the denominator is very far from converged, and the discontinuity factor is very far from one, which leads to a few possible observations: If one uses discontinuity factors early in the calculation, it could help speed up convergence since part of the reference solution is inserted. Unless we start with converged boundary fluxes, or a good guess, using condensed boundary fluxes for example, the discontinuity factor cannot be determined in a single transport sweep.
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
13
Fig. 3. Graphical representation of how to obtain the jump conditions in MOC and how to apply them.
4. Test case We implemented the proposed new method in a 2D MOC toy code. This 2D code was compared with OpenMOC (Boyd et al., 2014) to verify accuracy. The test-case is a square 2D pin-cell with a pitch of 1.26 cm and a pin radius of 0.4 cm unless stated otherwise. There is resonant absorption in the fuel, but resonance scattering is ignored in the initial development. The energy group width is 4 eV unless mentioned otherwise, and the resonance peak is located at 1000 eV in the center of the energy group. The cross section data is summarized in Table 1. The bare fuel pin model needed 128 azimuthal angles (Legendre quadrature), a track spacing of 0.001 cm and 3 polar angles (TY quadrature (Yamamoto et al., 2007)) to get results converged to a few pcms. For a lattice, the number of tracks and azimuthal angles needed is less as the tracks sample the pins multiple times when they cross the geometry. 5. Implementation and results At this point, the proposed method uses a discontinuity factor (DF) per angle and per region. This strategy has the same memory requirement as using angular cross sections, which would be preferred to using DFs. The numerator and denominator of these discontinuity factors can be averaged over groups of angles and then the same factor can be used for all angles in the group. This lowers significantly the memory cost. In Table 2, we can see the error on the fuel absorption reaction rates without equivalence (column 2) and with equivalence (column 3 to 9). On the left part of the table, we have applied discontinuity factors on the outgoing track angular fluxes only, and on the right side on the incoming and outgoing track angular fluxes. Then, in each of these two parts, we have allowed the discontinuity factors to depend either on: -azimuthal and polar angles and regions or -polar angles and regions only or -regions only (no angular dependence). When not depending on an angle set, we simply average the discontinuity factors over the angle set using the rele-
Table 1 Densities and non-resonant cross section data for the test case. Material
q(g/cc)
Nuclide
rsc: (b)
rabs: ðbÞ
UO2 H2O
10.4 0.7
U-238 O-16 H-1
11.29 3.888 20.47
0 0 0
vant quadrature. We can see that applying discontinuity factors on the incoming and outgoing fluxes yields the best results. When applying jump conditions depending on polar angles only, we do not lose accuracy compared with having azimuthally and polar angle dependent discontinuity factors. The discontinuity factors can have a polar and an azimuthal dependency, but a factor per surface independent of angle is enough to preserve partial currents. The polar dependence would be a unique feature compared to traditional equivalence methods because it accounts for axial effects. The reader should note that the errors mentioned are fractional errors on the fuel absorption rate for a single resonance in a single fuel pin. When scaling up to a full reactor problem, the maximum error coming from the flux separability approximation on U 238 absorption rates is closer to 2% and the mean error on eigenvalue closer to 400 pcm (Boyd, 2017). After the number of angles, the second factor of great memory costs is the number of regions. To decrease the memory requirements, we could make discontinuity factors a function of their neighborhood (Herrero et al., 1997), possibly with an artificial nodal network (Kozlowski et al., 2004), we could use clustering algorithms for discontinuity factors, like for the group crosssections (Boyd, 2017), or we could seek ways to propagate the factors inside a fuel pin from data gathered over a single pin region.
5.1. Parametric study At this point we want to show that this new method can be applied to a variety of cases derived from the resonant pin-cell problem. For that we perform a parametric study over pin radius, moderator density, presence of resonant scattering in the fuel and finally energy group width. Given the previous section’s results, we will only study polar angle dependent jump conditions applied on the incoming and outgoing track angular fluxes. So far, they are the best compromise between memory usage and accuracy of the reaction rates, with some potential for 3D equivalence due to the polar angle dependence.
5.1.1. Varying the fuel pin radius With increasing pin radius, the amount of self-shielding increases, so errors on the fuel absorption cross section should be larger. Table 3 shows results for varying pin radii. Without an equivalence method, the error first gets larger as expected without the correction, then reduces. This reduction is due to the compet-
14
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
Table 2 Fuel absorption reaction rate fractional errors (pcm) for different dependencies of the partial currents-based DFs. The error on the moderator flux is of the same order of magnitude. Discontinuity factors need to be applied on incoming and outgoing angular fluxes to preserve reaction rates. An angular dependence of DFs is not required to achieve conservation of reaction rates. DF used on:
No DF
Outgoing angular fluxes
DF dependencies:
No DF
Surface polar azimuthal
Surface polar
Surface
Incoming and outgoing angular fluxes Surface polar azimuthal
Surface polar
Surface
25 238 828 2206
6 23 47 29
3 15 32 38
9 47 115 153
1 1 1 0
0 0 0 0
0 0 0 0
rpeak ðbÞ 1e2 1e3 1e4 1e5
ing effect of the slowing down source in the fuel becoming more and more important as the fuel radius increases. For the discontinuity factors, first let’s observe that they are inducing much more correction on the pin-outgoing track angular fluxes than on the incoming and that the factors are usually above one. As previously stated, this is because there is insufficient selfshielding in the fuel when using the group cross-sections without their spatial and angular dependence. In order to reduce the group flux, the outgoing current has to be increased, which is why the factors are mostly above one. The discontinuity factors follow mostly the same trends as the error without equivalence. Then, we can see that the steepest and flattest polar angles exhibit opposite trends. The most horizontal polar angle discontinuity factors increases with pin radius, whereas the most vertical decreases with pin radius. This is driven by the source term in the fuel flux. The first/steepest polar angle is the one with the longest path in the fuel, so the contribution of the source is the largest for this angle. On the other hand, the length in the fuel along a track with the third polar angle is the shortest, so the contribution of the source compared to the contribution of the incoming flux is not as big. As mentioned above, the bigger the radius, the more the source in the fuel contributes to the fuel absorption reaction rate, so the lower the required correction. The correction introduced by the jump conditions is much larger for smaller pin radii. For larger pin radii, the distance ‘‘traveled” by the tracks in the moderator between exiting and entering gets smaller, meaning the effect of applying the outgoing jump condition is less dampened by the coolant before re-entering the fuel pin. These two effects could each introduce numerical instability. By numerical instability, we mean that the track angular fluxes could diverge because of consecutive multiplications by the discontinuity factors. As shown in this table, the method has reached convergence in all cases and reduced the error on reaction rates and group fluxes by orders of magnitude. 5.1.2. Varying the moderator density In order to show the applicability of our method for both boiling and pressurized water reactors, at zero power as well as full power, we vary the moderator density and the fuel pin radius simultaneously. Results shown in Table 4 indicate that the accuracy of our
method is not affected. As expected, the higher the water density, the stronger the source in the water and the higher the error in the absorption rate due to the error in self-shielding when using the group-cross section without equivalence factors. 5.1.3. Introducing resonant scattering in the fuel So far we have only modeled resonant absorption in the fuel. A more realistic model would include resonant scattering. This could create issues with the equivalence, as now the error comes both from the condensation of the absorption and scattering crosssections. In order to capture the up-scatter and down-scatter rates of neutrons properly it is important to get the scattering rates right. Our method preserves group fluxes, so it preserves all reaction rates indiscriminately. To further support this claim, we present in Table 5 results with and without resonant scattering in the fuel. In order to introduce scattering without changing the absorption cross section energy profile, we choose to adapt the resonance widths of scattering and absorption. The total width is kept at 0.025 eV, while the absorption width is 0.015 eV and the scattering width 0.01 eV. The resonant scattering energy dependence is then exactly the same as in Fig. 1a. We first observe that the introduction of scattering diminishes the absorption group cross section. In the narrow resonance model, when neutrons scatter, they scatter out of the resonance and can’t be absorbed anymore. Even though there is less absorption going on, the error when not using equivalence methods increases in the presence of scattering. Resonant scattering in the narrow resonance model acts very similarly as absorption, as only potential scattering contributes to the source term. So the self-shielding effectively increases with the presence of scattering, and the error on absorption rates gets larger. Finally, as the self shielding increases, the equivalence factors are getting larger to correct the increasing error on the group fluxes and reaction rates. 5.1.4. Varying energy group widths Because the effect of resonant absorption on the absorption group cross section strongly depends on the energy group width, we need to make sure that our method is applicable for the wide range of energy group widths that is used in current deterministic transport solvers. In Table 6, the method is tested for groups of
Table 3 Fractional errors on fuel absorption reaction rate and polar discontinuity factors for various radii of a fuel pin with a 104 b peak SLBW absorption-only resonance in the fuel. DFs on the incoming angular fluxes only perform minor corrections, and DFs on the outgoing angular fluxes are larger for cases with larger errors without equivalence. Polar angle 1 is the steepest in the quadrature. Pin radius (cm)
Error on fuel absorption rate (pcm)
Polar incoming DF
Polar outgoing DF
No DF
DF
1
2
3
1
2
3
0.2 0.3 0.4 0.5 0.6
859 910 829 621 261
<0.1 <0.1 <0.1 <0.1 <0.1
1.00001 1.00002 1.00009 1.00027 1.00067
1.00011 1.00021 1.00032 1.00031 1.00051
1.00024 1.00019 1.00010 0.99995 1.00017
1.14989 1.13096 1.10441 1.07431 1.04102
1.02155 1.02266 1.02103 1.01684 1.01028
0.99671 0.99827 0.99999 1.00132 1.00167
15
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
Table 4 Fractional errors on fuel absorption reaction rate for various radii of a fuel pin at various water densities with a 104 b peak SLBW absorption-only resonance in the fuel. The pin-topin pitch is still 1.26 cm. DFs allow for exact preservation of reaction rates in all cases, representative of LWR operating conditions. Fractional error on fuel absorption rate (pcm) 0.2 g/cc Pin rad.
0.7 g/cc
1 g/cc
No DF
DF
No DF
DF
No DF
DF
643 580 401 209 50
<0.1 <0.1 <0.1 <0.1 <0.1
859 910 829 621 261
<0.1 <0.1 <0.1 <0.1 <0.1
896 979 936 754 363
<0.1 <0.1 <0.1 <0.1 <0.1
qmod 0.2 0.3 0.4 0.5 0.6
Table 5 Fractional errors on fuel absorption reaction rate (in pcm), absorption group cross sections and polar outgoing discontinuity factors with resonant absorption and without/with resonant scattering. The introduction of scattering shields the absorption but also increases the self-shielding error. Discontinuity factors increase to provide additional correction and allow for exact preservation of reaction rates. Without resonant scattering
With resonant scattering
rpeak
rabs;g
Fuel abs. error
Polar outgoing DF
(b)
(b)
No DF
DF
Angle 1
2
1e1 1e2 1e3 1e4 1e5
.09 .63 2.4 8.3 30
0.4 25 238 828 2206
<0.1 <0.1 <0.1 <0.1 <0.1
1.00029 1.00703 1.03362 1.10441 1.23649
1.00002 1.00084 1.00567 1.02103 1.08287
rabs;g
Fuel abs. error
Polar outgoing DF
3
(b)
No DF
DF
Angle 1
2
3
0.99999 0.99983 0.99954 0.99999 1.01376
.09 .55 2.0 6.8 25
3.2 60 357 1168 3059
<0.1 <0.1 <0.1 <0.1 0.9
1.00205 1.01590 1.06050 1.19302 1.49169
1.00014 1.00188 1.00881 1.03188 1.13024
0.99995 0.99964 0.99913 0.99917 1.01872
Table 6 Fractional errors on fuel absorption reaction rate and polar discontinuity factors for various energy group-widths (eV) with a 104 barns peak SLBW absorption-only resonance in the fuel. The smaller the energy group width, the higher the error on the absorption reaction rate, and the larger the discontinuity factors. Group width (eV)
20 10 4 2 1 0.5
Fuel abs. error (pcm)
Polar incoming DF
No DF
DF
1
2
3
1
2
3
174 345 829 1535 2504 2737
<0.1 <0.1 <0.1 <0.1 <0.1 <0.1
1.00002 1.00005 1.00009 1.00010 1.00005 0.99999
1.00008 1.00015 1.00032 1.00052 1.00070 1.00055
1.00003 1.00005 1.00010 1.00013 1.00017 1.00048
1.02241 1.04416 1.10441 1.18303 1.24798 1.17955
1.00381 1.00781 1.02103 1.04687 1.10906 1.21131
0.99972 0.99957 0.99999 1.00373 1.02396 1.11074
widths ranging from 20 to 0.5 eV, centered on a single SLBW resonance. We can see that the condensation error without an equivalence method gets larger as the group width decreases. Simultaneously, the discontinuity factors increase mostly monotonously to correct the partial currents.
6. Conclusion and future work Cross sections for multi-group transport calculations are traditionally generated using the scalar flux, thus neglecting the angular dependence of group cross sections. This approximation leads to a large systematic bias on pin power distributions in light water reactor systems, and an equivalence method has to be used to preserve reaction rates. This paper introduces a new state-of-the-art equivalence method that incorporates an extended level of transport solution conservation; namely direct polar angle current preservation. It is based on jump conditions applied on the incoming and outgoing track angular fluxes at the fuel pin surface. These jump conditions are chosen to be polar angle dependent to enforce polar angle partial currents preservation. This feature is expected to be important for 3D MOC transport calculations. The method is tested on a 2D fixed source pin-cell problem by condensing fine energy resonance fluxes and resonant cross sec-
Polar outgoing DF
tions in a parametric study on pin radius, moderator density, presence of resonant scattering and energy group width. The parameter ranges were chosen to encompass LWR operating conditions, and common group widths in deterministic multi-group transport. The method demonstrates precise preservation of all reaction rates in all cases. This equivalence method has already been extended to the case of spatial representations with multiple fuel rings commonly used in fuel depletion and fuel temperature distribution modeling for LWRs and also showed exact preservation. While this article assumed reference partial currents were known in order to derive equivalence parameters, future work will focus on alternative ways to obtain the jump conditions. This work will also soon be extended to the case of Monte Carlo-generated group cross sections so that one can obtain accurate equivalence parameters directly from partially converged, full-core, Monte Carlo calculations; thereby avoiding all approximations inherent in traditional multi-group cross section generation methods.
Acknowledgements The research work presented in this paper was supported by the US Department of Energy Nuclear Energy University Program contract DE-NE0008578.
16
G. Giudicelli et al. / Annals of Nuclear Energy 112 (2018) 9–16
References Bell, G.I., Hansen, G.E., Sandmeier, H.A., 1967. Multitable treatments of anisotropic scattering in SN multigroup transport calculations. Nucl. Sci. Eng. 379, 28. Boyd, W., 2017. Reactor Agnostic Multi-Group Cross Section Generation for FineMesh Deterministic Neutron Transport Simulations (Ph.D. thesis). Massachusetts Institute of Technology. Boyd, W., Gibson, N., Forget, B., Smith, K., 2017. An analysis of condensation errors in multi-group cross section generation for fine mesh neutron transport calculations. Submitted to Annals of Nuclear Energy. Boyd, W., Shaner, S., Li, L., Forget, B., Smith, K., 2014. The openmoc method of characteristics neutral particle transport code. Ann. Nucl. Energy 68, 43. Gibson, N., Smith, K., Forget, B., 2015. Simple benchmark for evaluating selfshielding models. In: Proceedings of M&C 2015, Nashville, Tennessee. Grassi, G., 2007. A nonlinear space-angle multigrid acceleration for the method of characteristics in unstructured meshes. Nucl. Sci. Eng. 155 (2), 208–222. Gunow, G., Shaner, S., Boyd, W., Forget, B., Smith, K., 2017. Accuracy and performance of 3D MOC for full-core PWR problems. In: Proceedings of M&C2017, Jeju Korea. Hebert, A., Benoist, P., 1991. A consistent technique for the global homogenization of a pressurized water reactor assembly. Nucl. Sci. Eng. 109 (4), 360–372. Herrero, J., Garcia-Herranz, N., Ahnert, C., 1997. Neighborhood corrected interface discontinuity factors for multi-group pin-by-pin diffusion calculations for LWR. Ann. Nucl. Energy 24 (6), 477–488.
Kavenoky, A., 1978. The SPH homogenization method. In: IAEA Proceedings of Specialists Meeting on Homogenization Method in Reactor Physics, Lugano, Switzerland. Koebke, K., 1978. A new approach to homogenization and group condensation. In: IAEA Proceedings of Specialists Meeting on Homogenization Method in Reactor Physics, Lugano, Switzerland. Kozlowski, T., Lee, D., Downar, T.J., 2004. The use of artificial neural network for online prediction of pin-cell discontinuity factors in PARCS. In: PHYSOR 2004 Proceedings. Kozlowski, T., Xu, Y., Downar, T.J., Lee, D., 2011. Cell homogenization method for pin-by-pin neutron transport calculations. Nucl. Sci. Eng. 169, 1–18. Miao, J., Forget, B., Smith, K., 2016. Analysis of correlations and their impact on convergence rates in Monte Carlo eigenvalue simulations. Ann. Nucl. Energy 92, 81–95. Sanchez, R., 2009. Assembly homogenization techniques for core calculations. Prog. Nucl. Energy 51, 14–31. Smith, K., 1984. Nodal method storage reduction by non-linear iteration. Trans. ANS 44, 265. Yamamoto, A., Tabuchi, M., Sugimura, N., Ushio, T., Mori, M., 2007. Derivation of optimal polar angle quadrature set for the method of characteristics based on approximation error of the Bickley functions. J. Nucl. Sci. Technol. 44 (2), 129– 136.