- Email: [email protected]
COMET neutronics solutions to the prismatic VHTR benchmark problem
Nuclear Engineering and Design 358 (2020) 110395
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
COMET neutronics solutions to the prismatic VHTR benchmark problem Dingkang Zhang, Farzad Rahnema
⁎
T
Georgia Institute of Technology, 770 State Street NW, Atlanta, GA 30332-0745, USA
ARTICLE INFO
ABSTRACT
Keywords: Hybrid transport methods Incident flux response expansion Benchmark calculations Very high temperature reactors
Recently, the hybrid stochastic deterministic coarse mesh radiation transport method COMET for whole core criticality calculations was extended to treat the energy variable continuously consistent with its treatment of the rest of the phase space variables. In this paper, a surface-dependent asymptotic spectrum expansion scheme was used by the continuous-energy COMET to obtain the transport solutions to a stylized Very High Temperature Reactor (VHTR) benchmark problem in which all levels of heterogeneity that are of reactor physics importance are preserved. Because of its complex geometry and material distribution, this benchmark problem is challenging to solve with deterministic neutronic methods, even those based on transport theory. The eigenvalue and the fuel pin fission density distribution for all control rods fully inserted and fully withdrawn configurations were compared to the corresponding values calculated by the continuous-energy Monte Carlo code MCNP. The agreement between the two codes was found to be excellent. The eigenvalue difference was 66 pcm for the uncontrolled configuration and 69 pcm for the controlled configuration. The corresponding average pin fission density differences were 0.36% and 0.53% for the uncontrolled and controlled configurations, respectively. It was also found that COMET is 4–5 orders of magnitude faster than MCNP. This indicates that COMET has continuous-energy Monte Carlo like accuracy but with formidable computational speed for solving problems with multiple levels of heterogeneity such as the VHTR problem.
1. Introduction The Very High Temperature Reactor (VHTR) core design is very challenging to solve by any neutronic methods. This is mainly due to the multiple levels of heterogeneity resulting from the use of TRISO fuel particles embedded in a graphite matrix which includes burnable absorber particles in the corners of the fuel assemblies as well as control rods in the graphite reflector (non-fissionable) blocks. Stochastic transport methods regarded as the gold standard require significant computational resources to generate detailed solutions (e.g., fuel pin power distribution) that are statistically converged. Performance of deterministic transport methods heavily relies on the accuracy of the multigroup cross sections which are precomputed by local (e.g., single block) transport calculations. However, this accuracy suffers from the approximations made during the cross section processing. These include spatial homogenization (e.g., of the fuel and burnable absorber particles), control rod cross section generation in non-fissionable regions (reflector blocks), and approximating the core environment effect that is significant because the local spectrum is driven by the neighboring fuel and reflector blocks (i.e., the local boundary conditions are not known a priori). Low order transport methods such as diffusion
⁎
theory suffer from the same issues on top of the need for a broader homogenization (e.g., in an entire fuel block and reflector blocks with and without control rods) and the use of the usual approximations made to derive the theory. Therefore, it is prudent to assess the accuracy of the neutronics methods that are used for VHTR calculations. In this paper, we evaluate the accuracy of COMET for VHTR calculations using the prismatic VHTR benchmark problem (Connolly et al., 2015). COMET is a hybrid stochastic/deterministic transport method based on the incident flux response expansion theory (Mosher and Rahnema, 2006). The numerical implementation of the method (Zhang and Rahnema, 2012) consists of using response expansion coefficients (RC) generated by local calculations followed by core calculations to obtain the desired solution. A stochastic method is used to pre-compute a library of the incident flux RCs for all the unique coarse meshes (such as fuel and reflector blocks) in the core. Using this library, deterministic sweeps are performed to converge on the outgoing/incoming fluxes on coarse mesh interfaces in the core. The desired whole core solution is then generated by a superposition of the converged incident responses and the precomputed RCs. In the original method (Zhang and Rahnema, 2012), the phase space was treated continuously except for the energy variable. Extensive benchmarking of COMET
Corresponding author. E-mail address: [email protected] (F. Rahnema).
https://doi.org/10.1016/j.nucengdes.2019.110395 Received 11 July 2019; Received in revised form 15 October 2019; Accepted 21 October 2019 0029-5493/ © 2019 Published by Elsevier B.V.
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
(Zhang and Rahnema, 2017; Ulmer et al., 2016; Connolly and Rahnema, 2013) against the Monte Carlo method in various 3D problems typical of different reactor types has confirmed its high fidelity and formidable computational efficiency. Recently, COMET was extended (Rahnema and Zhang, 2018) to also treat the energy variable continuously. Like stochastic methods, this extension allows COMET to fully model the material heterogeneity down to the fuel and burnable absorber particle levels. As a result, in addition to any reactor type, it is demonstrated in this paper that COMET is an ideal tool for modeling the VHTRs because of its continuous-energy stochastic-like fidelity while retaining the formidable computational efficiency. The paper is organized as follows. For the sake of completeness and the reader’s convenience, the continuous energy COMET method is briefly reviewed in Section 2 and the VHTR benchmark problem is summarized in Section 3. The numerical results are presented and discussed in Section 4 and the paper is concluded in Section 5.
Rism ( r ,
Jis+, m =
1 F (r , k i
, E ) for r
Vi ,
,
, E) =
i
(r ,
Rismsm =
,
Jis±, m W (E )
, E) =
m( r
,
, E ) for r
m
m( r
(3)
In the above equation, represent the outgoing or incoming angular flux on surface s of coarse mesh i, denoted as Vis , m represents a set of orthogonal expansion functions, W (E ) is an asymptotic spectrum (Rahnema and Zhang, 2018), and Jis± , m are the corresponding expansion coefficients. Based on the incident flux response expansion theory, the flux distribution in each coarse mesh can be computed as the following superposition. ± is
i(r
,
Jis , m Rism ( r ,
, E) =
, E )for r
s, m
Vi .
, E) =
1 FRism ( r , k
, E)
s'm'
, m'
'
R m' m , is s
(7)
Vis
dr
nis+·
>0
d
dE (nis+· ) Rism ( r ,
, E)
m( r
,
, E)
,
, E ) = Pl (x ) Pn (y ) Up (µ ) Pq ( ) Bg,0 (E ),
(9)
3. Stylized prismatic VHTR benchmark problem A brief description of the stylized VHTR benchmark problem is presented in this section. Readers are referred to Connolly et al. (2015) for the detailed description. As shown in Fig. 1, the VHTR benchmark problem consists of a graphite inner reflector, three-rings of fuel blocks (yellow regions in Fig. 1), and an outer graphite reflector. The green regions in the reflectors represent the locations of the control channels, while the corresponding regions in the fuel blocks denote the reserved shutdown channels. In the axial direction, the core consists of a two-layer bottom reflector, a ten-layer active core, and a two-layer bottom reflector. The height of each layer is 79.3 cm. Vacuum boundary condition is imposed on all external surfaces. The following two core configurations are modeled in this paper. We note that a critical core configuration was not provided for the continuous energy case.
(4)
Here, Rism are the surface-to-volume response expansion functions that can be obtained by solving the following local fixed-source problem:
HRism ( r ,
is'
where the functions P and U are the Legendre polynomials and the Chebyshev polynomials of the second kind, respectively, and B represents the B-splines. That is, Legendre polynomials are used to expand the spatial dependence and a tensor product of Chebyshev polynomials of the second kind and Legendre polynomials are used to expand the polar and azimuthal angular distribution, respectively. It can be seen from Eqs. (3) and (9) that the energy distribution is expanded in terms of the products of the asymptotic spectrum W(E) and the 0th order B-splines with coarse energy bins. The numerical implementation of the COMET method consists of the following steps. In the first step, a continuous energy stochastic method is used to solve the local fixed-source problems defined by Eqs. (5) and (6) to obtain the library of the incident response coefficients. Because of the geometric robustness and the continuous energy treatment even the heterogeneity within the fuel and burnable absorber particles can be fully modeled. This detailed treatment of material heterogeneity is performed through the stochastic precomputation of the response coefficient (function) library. The precomputation of the library is similar to the industry practice of precomputation of the multigroup cross section library. In the second step, an iterative deterministic sweeping solver (Zhang and Rahnema, 2012) is used to converge the core solution via Eq. (7). In the last step, Eq. (4) is used to construct the flux (or the pin power) distribution in each coarse mesh by superposition.
(2)
Vis
J
In theory, any complete function set can be used in Eq. (3) to expand the outgoing/incoming angular fluxes on coarse mesh surfaces as long as they satisfy the orthogonal conditions. However, extensive benchmarks (Zhang and Rahnema, 2017; Ulmer et al., 2016; Connolly and Rahnema, 2013) have shown the following expansion functions are a good choice since they can sufficiently represent the angular flux phase distributions at a relatively low order and consequently achieve high computational efficiency while maintaining high accuracy.
where H and F are the standard transport and fission operators, k stands for the core eigenvalue (keff), i represents the angular flux in coarse mesh i at location r or (x , y , z ) in the direction of the unit vector or (µ, ) and with energy E , i is the flux from the neighboring meshes incident on the boundary of mesh i ( Vi ), and ni+ denotes the outward normal to Vi . The variables µ and are the cosine of the polar angle and the azimuthal angle, respectively. Since the incident flux i on each surface is not known a priori, it is assumed it can be expanded in terms of a set of orthogonal expansion functions as below. ± is ( r
(6)
Vis & · nis+ < 0
(8)
(1)
Vi & ·ni+ < 0,
, E )for r
, E )for r
where R m' m is the surface-to-surface response coefficients that can be is s computed by projecting the sourface-to-volume response functions m Ris ( r , , E ) onto the expansion basis.
with the boundary condition i(r
,
'
The detailed description of the COMET method and the continuous energy treatment can be found in Mosher and Rahnema (2006), Zhang and Rahnema (2012) and Rahnema and Zhang (2018). A short description is presented in this section for the sake of completeness and the readers convenience. To solve the Boltzmann transport equation in a reactor core, the spatial domain of the problem, V , is divided into a number of nonoverlapping coarse meshes {Vi i = 1, 2, , N } , and the phase space distribution of the flux in each coarse mesh can be obtained by solving the following set of local problems.
, E) =
m( r
The partial current expansion moments Jis+, m can be iteratively determined by projecting Eq. (2) onto the expansion basis as below.
2. Continuous energy COMET method
H i(r ,
, E ) = W (E )
(5)
with the boundary condition 2
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Fig. 1. Radial configuration of the VHTR benchmark problem; Taken from Connolly et al. (2015).
Fig. 3. Cut through, fuel pin; Taken from Connolly et al. (2015).
• All rods out (ARO): all control rods are withdrawn. • All rods in (ARI): all control rods in the inner and outer reflectors are
(assemblies) consist of small spheres (with a radius of 0.02125 cm) of uranium oxycarbide enriched to 15% in 235U by weight surrounded by four layers of materials: carbon, pyrolitic carbon, silicon carbide, and pyrolitic carbon. The packing fraction within the fuel pins is 35%. These fuel particles are modeled as a body centered cubic (BCC) lattice with the side length of 0.09664222 cm. it is noted that the choice of the BCC leads to unrealistic case of cutting particles/spheres as seen in Fig. 3. This of course is inconsequential from a benchmarking point of view.
fully inserted.
3.1. Fuel blocks The geometric configuration of the fuel blocks (with or without the reserve shutdown channel) is illustrated in Fig. 2. The apothem of each fuel block is 18 cm. Each fuel block consists of a hexagonal graphite structure (grey area in Fig. 2), 210 fuel pins (yellow regions), 108 coolant channels (white regions), 6 burnable absorber pins (green regions). Fuel blocks with the reserve shutdown channel have a similar geometric configuration as the standard fuel blocks except they have a 9.525-cm diameter hole located 9.756-cm away from the center. As a result, they have 186 fuel pins, 95 coolant channels, and 6 burnable absorber pins.
4. Numerical results The benchmark problem was modeled by both COMET and MCNP (X-5 Monte Carlo, 2005). The eigenvalues and the fission density distributions predicted by the two codes were then compared. The following four statistical quantities as defined in Zhang and Rahnema (2012) were used to summarize the differences between the pin fission density values (distributions) predicted by the two codes: the average (AVG), root mean square (RMS), mean relative/flux weighted (MRE), maximum (MAX) of the relative differences in the pin fission densities.
3.2. Fuel pins Fuel pins (see Fig. 3) imbedded in the graphite fuel blocks
Fig. 2. Configuration of fuel blocks and fuel blocks with the reserve shutdown channel; Taken from Connolly et al. (2015).
3
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Fig. 4. MCNP computation model.
4.1. MCNP calculation Because the core configurations have 1/12th symmetry in the radial direction and 1/2 symmetry in the axial direction, only 1/24th of the core was modeled in MCNP as shown in Fig. 4. This significantly reduced the MCNP computation time to maintain the same uncertainty as in the full core case. More importantly, this also avoided the MCNP fission source convergence problem that typically is a very challenging issue for large eigenvalue problems as fully discussed/analyzed in Zhang and Rahnema (2017). MCNP calculations were performed on an 84-CPU cluster. The fission source distribution was first converged by following 50 million particle histories, and then 25 million inactive and 500 million (4000 active cycles with 125,000 histories in each cycle) active particle histories were used to perform the actual calculations for each core configuration.
Fig. 5. Color set consisting of fuel and graphite blocks for generating asymptotic spectra.
surfaces. In the second approach, the spectrum on surface 5 from the color set calculation was used as the asymptotic function W (E ) on all fuel-fuel surfaces, the average of the incoming/outgoing spectra on surfaces 3 and 4 were used as the asymptotic functions on all fuel-reflector/reflector-fuel surfaces, and the average of the spectra on surfaces 1 and 2 was used as the asymptotic function on all reflector-reflector surfaces. The stochastic continuous energy response function generator (modified MCNP) was then used to pre-compute the response coefficients for all unique coarse meshes. The expansion order was set to {4,4,2,2} which represents the 4th order expansion on the two spatial variables on coarse mesh boundaries (surfaces) and the 2nd order expansion on the two angular (polar and azimuthal) variables. The energy distribution deviation away from W(E) was expanded by the 0th order B-splines with eight coarse energy bins with bin boundaries of 1.4572E−7, 6.2506E−7, 3.9279E−6, 1.3007E−4, 9.1188E−3, 8.2085E−1, 2.2313, and 20 MeV. It should be noted that the energy bin structure was originally used for the Pressurized Water Reactor benchmark problems (Rahnema et al., 2013) and may not be optimal for the VHTR problems. Regardless, it is demonstrated that the method works very well despite a choice that is not necessarily optimized or ideal. Of course, optimization of the energy bin structure is expected to further improve the COMET accuracy. However, this study is beyond the scope of this work. The response function library was precomputed with 15 million particle histories in each calculation that took about 3.2 h for a fuel block and 6 min for the other coarse meshes on a single CPU. The total computation time for the response function generation was about 873 h, or equivalently 10.4 h on the 84-CPU cluster.
4.2. COMET calculations Nine unique coarse meshes were used by COMET to model the VHTR benchmark problem. They include: two fuel blocks (without and with the reserve shutdown channel), two cooled reflector blocks (without and with the control rod hole), 3 hexagonal graphite reflector blocks (without the control rod hole, with an empty control rod hole, and with the control rod inserted), and two incomplete graphite blocks to match the core external boundary as shown in Fig. 1. The dimension and properties of each coarse mesh are given in Table 1. It can be seen from Section 2 that the choice of the asymptotic spectrum W (E ) is not unique. However, if W (E ) is very close to the actual spectrum, it can significantly lower the energy expansion order. In this work, the following two choices were tested: 1) a single asymptotic spectrum for all surfaces, and 2) surface-dependent asymptotic spectra based on its location in the core. A color set shown in Fig. 5 consisting of three graphite reflector blocks (coarse mesh 5 in Table 1) and three fuel blocks (coarse mesh 1 in Table 1) with the specular boundary condition was first used to generate the asymptotic spectra. In the first approach, the energy distribution on surface 5 (Fig. 5) was used as the asymptotic spectrum on all coarse mesh Table 1 Dimension and property of each unique coarse mesh. Index
Coarse mesh
Dimension: Apothem × Height (cm)
1 2 3 4 5 6 7 8 9
Fuel block Fuel block with reserve shutdown channel Cooled reflector Cooled reflector with reserve shutdown channel Reflector block Reflector block with control rod hole Reflector block with control rod inserted Side reflector Corner reflector
18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 1/2 of 18 × 49.3 hexagon 1/3 of 18 × 49.3 hexagon
4
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Table 2 Comparison of the core eigenvalues and computation times. Core
ARO
ARI
Method
COMET**
MCNP*
Eigenvalue Difference (pcm) Run Time (hours)
1.06617 ± 0.00004
Eigenvalue Difference (pcm) Run Time (hours)
0.91534 ± 0.00004
49.2
49.5
Single Asymptotic Spectrum
Surface-dependent Asymptotic Spectrum
1.06464 ± 0.00005 153 ± 6 0.89
1.06683 ± 0.00005 66 ± 6 0.90
0.91674 ± 0.00005 140 ± 6 0.85
0.91603 ± 0.00005 69 ± 6 0.85
* MCNP calculations modeled only 1/24th of the core on an 84-CPU cluster. It is estimated that a single full core calculation on this cluster would take about 50 days. ** COMET calculations modeled the full core on a single CPU of the cluster.
The response function library was used to perform core calculations for the two core configurations. To demonstrate the COMET computational efficiency, the full core was modeled in COMET. The core eigenvalue and the fuel pin fission density distributions computed by COMET are compared to the corresponding MCNP results in Section 4.3.
plotted in Figs. 6 and 7. Because of the axial symmetry only the top five axial layers are presented. For the reader’s interest, the pin fission density distributions predicted by COMET for the ARO and ARI cases are plotted in the appendix (Figs. A1 and A2), respectively. It can be seen from Figs. 6 and 7 that the maximum differences in the pin fission densities occurs in the peripheral fuel pins in the top layer next to the inner and outer reflectors when the single asymptotic spectrum was used. This is because a single asymptotic spectrum W (E ) with coarse energy bins is not sufficient for representing all surface spectra such as those on the fuel-fuel, fuel-reflector, reflector-reflector, reflector-vacuum, etc. Such an approximation leads to an underestimation of the thermal neutron flux on interfaces between the graphite blocks and the fuel blocks. As seen from Figs. 6b and 7b, the agreement between COMET and MCNP is excellent when the surface dependent asymptotic spectra were used on each surface without compromising the formidable computational efficiency of COMET. Note that the response function pre-computation time increased by a factor 2.
4.3. Comparison of calculated results The core eigenvalues predicted by COMET and MCNP and the corresponding computation times are presented in Table 2. It can be seen that the differences between the eigenvalues predicted by COMET with the single asymptotic spectrum and MCNP are 153 pcm and 140 pcm for the ARO and ARI cases, respectively. The corresponding differences for the surface-dependent spectra are 66 pcm and 69 pcm. The COMET full core calculations took only about 0.9 h on a single core, while the MCNP calculation for both of the 1/24th-core configurations took about 49 h each on the 84-CPU cluster as seen in Table 2. The comparison of the pin fission density distribution predicted by COMET and MCNP are statistically summarized in Table 3. The four statistical parameters are defined as follows: the average (AVG), the mean relative (MRE), the root mean square (RMS) and the maximum (MAX) relative differences. The average and maximum relative differences between COMET with the single asymptotic spectrum and MCNP are 1.04% and 3.72% for the ARO case and 0.98% and 4.60% for the ARI case, respectively. When the surface-dependent asymptotic spectra are used then the agreement between COMET and MCNP is improved significantly as expected. The corresponding differences are 0.36% and 1.84% for the ARO case and 0.53% and 2.13% for the ARI case. The difference between the COMET and MCNP pin fission density distributions are
5. Summary and conclusions In this paper, the performance of the recently developed continuous energy COMET method was evaluated for modeling VHTR problems. This was done by comparing the COMET solution (core eigenvalue and detailed fuel pin fission density distribution) to the corresponding MCNP solution in the stylized prismatic VHTR benchmark problem with control rods fully inserted and withdrawn. COMET treats the energy variable continuously by decomposing the angular flux incident on the coarse mesh surfaces as a product of a fast-varying function in energy, called the asymptotic spectrum, and a slowly varying function in energy. It was found that the use of a (color) set of asymptotic spectra (namely, fuel-fuel, fuel-reflector, reflector-fuel, and reflector-reflector) leads to excellent results. The differences between the eigenvalues calculated by COMET and MCNP were 66 pcm for the ARO case and 69 pcm for the ARI case. The average and maximum relative difference in the pin fission density distributions were 0.36% and 1.84% for the ARO case and 0.53% and 2.13% for the ARI case, respectively. COMET was also found to be significantly faster than MCNP. Each COMET full core calculation took about 0.9 h on a single CPU, while MCNP calculations for 1/24th of the core on an 84-CPU cluster took about 49 h each. It should be noted that if the full-core problems were to be modeled then at least 24 times more particle histories would have been needed to maintain the same statistical uncertainty. However, obtaining a fully symmetric (converged) solution will require averaging solutions from multiple (50–100) independent runs (Zhang and Rahnema, 2017). Conservatively, assuming a scaling factor of 1, this leads to about 6.8 years – 13.6 years of MCNP computation time on the 84-CPU cluster. This indicates that COMET is 4–5 orders of magnitude faster than MCNP. It can be concluded that COMET is an ideal tool for
Table 3 The statistics of the differences between the Pin fission density distribution predicted by COMET* and MCNP**. Core
ARO
ARI
Asymptotic Spectrum
Single
Surface-dependent
Single
Surface-dependent
AVG MRE RMS MAX
1.04% 1.09% 1.23% 3.72%
0.36% 0.34% 0.46% 1.84%
0.98% 0.95% 1.20% 4.60%
0.53% 0.52% 0.64% 2.13%
* The COMET pin fission density uncertainties are in the range of 0.042% – 0.064%. ** The MCNP pin fission density uncertainties are in the range of 0.16% – 0.30%. 5
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Fig. 6. Relative difference between the pin fission density distributions in the top five axial layers predicted by COMET and MCNP for the ARO case.
Fig. 7. Relative difference between the pin fission density distribution in the top five axial layers predicted by COMET and MCNP for the ARI case.
modeling the VHTRs because of its continuous-energy stochastic-like fidelity but with a formidable computational speed.
Acknowledgements The second author owns equity in a company that has licensed the COMET technologies from Georgia Tech. This study which is a demonstration of COMET could affect his personal financial status. The terms of this arrangement have been reviewed and approved by Georgia Tech in accordance with its conflict of interest policies.
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 6
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Appendix A. . Comet pin fission density distributions The normalized pin fission density distributions in the top five axial layers predicted by COMET for the ARO and ARI cases are plotted in Figs. A1 and A2, respectively.
Fig. A1. Normalized pin fission density distribution in the top five layers predicted by COMET with surface-dependent asymptotic spectra for the ARO case.
Fig. A2. Normalized pin fission density distribution in the top five layers predicted by COMET with surface-dependent asymptotic spectra for the ARI case.
7
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
References
reactor benchmark problem with UO2 and MOX fuel. Nucl. Technol. 184, 1–28. Ulmer, R., Rahnema, F., Connolly, K., 2016. A neutronic benchmark specification and COMET solution for the Advanced Burner Test Reactor. Ann. Nucl. Energy 87, 76–106. X-5 Monte Carlo Team, 2005. MCNP—A General Monte Carlo N-Particle Transport Code, Version 5. Los Alamos National Laboratory. Zhang, D., Rahnema, F., 2012. An efficient hybrid stochastic/deterministic coarse mesh neutron transport method. Ann. Nucl. Energy 41, 1–11. Zhang, D., Rahnema, F., 2017. Monte Carlo solution convergence analysis via COMET in a stylized integral inherently safe LWR benchmark problem. Ann. Nucl. Energy 100, 3–11.
Connolly, K., Rahnema, F., 2013. A heterogeneous coarse mesh radiation transport method for neutronic analysis of prismatic reactors. Ann. Nucl. Energy 56, 87–101. Connolly, K., Rahnema, F., Tsvetkov, P., 2015. Prismatic VHTR neutronic benchmark problems. Nucl. Eng. Des. 285, 207–240. Mosher, S., Rahnema, F., 2006. The incident flux response expansion method for heterogeneous coarse mesh transport problems. Transp. Theory Stat. Phys. 34, 1–26. Rahnema, F., Zhang, D., 2018. Continuous energy coarse mesh transport (COMET) method. Ann. Nuclear Energy 115. Rahnema, F., Hon, R., Douglass, S., 2013. A stylized three-dimensional pressurized water
8
Contents lists available at ScienceDirect
Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes
COMET neutronics solutions to the prismatic VHTR benchmark problem Dingkang Zhang, Farzad Rahnema
⁎
T
Georgia Institute of Technology, 770 State Street NW, Atlanta, GA 30332-0745, USA
ARTICLE INFO
ABSTRACT
Keywords: Hybrid transport methods Incident flux response expansion Benchmark calculations Very high temperature reactors
Recently, the hybrid stochastic deterministic coarse mesh radiation transport method COMET for whole core criticality calculations was extended to treat the energy variable continuously consistent with its treatment of the rest of the phase space variables. In this paper, a surface-dependent asymptotic spectrum expansion scheme was used by the continuous-energy COMET to obtain the transport solutions to a stylized Very High Temperature Reactor (VHTR) benchmark problem in which all levels of heterogeneity that are of reactor physics importance are preserved. Because of its complex geometry and material distribution, this benchmark problem is challenging to solve with deterministic neutronic methods, even those based on transport theory. The eigenvalue and the fuel pin fission density distribution for all control rods fully inserted and fully withdrawn configurations were compared to the corresponding values calculated by the continuous-energy Monte Carlo code MCNP. The agreement between the two codes was found to be excellent. The eigenvalue difference was 66 pcm for the uncontrolled configuration and 69 pcm for the controlled configuration. The corresponding average pin fission density differences were 0.36% and 0.53% for the uncontrolled and controlled configurations, respectively. It was also found that COMET is 4–5 orders of magnitude faster than MCNP. This indicates that COMET has continuous-energy Monte Carlo like accuracy but with formidable computational speed for solving problems with multiple levels of heterogeneity such as the VHTR problem.
1. Introduction The Very High Temperature Reactor (VHTR) core design is very challenging to solve by any neutronic methods. This is mainly due to the multiple levels of heterogeneity resulting from the use of TRISO fuel particles embedded in a graphite matrix which includes burnable absorber particles in the corners of the fuel assemblies as well as control rods in the graphite reflector (non-fissionable) blocks. Stochastic transport methods regarded as the gold standard require significant computational resources to generate detailed solutions (e.g., fuel pin power distribution) that are statistically converged. Performance of deterministic transport methods heavily relies on the accuracy of the multigroup cross sections which are precomputed by local (e.g., single block) transport calculations. However, this accuracy suffers from the approximations made during the cross section processing. These include spatial homogenization (e.g., of the fuel and burnable absorber particles), control rod cross section generation in non-fissionable regions (reflector blocks), and approximating the core environment effect that is significant because the local spectrum is driven by the neighboring fuel and reflector blocks (i.e., the local boundary conditions are not known a priori). Low order transport methods such as diffusion
⁎
theory suffer from the same issues on top of the need for a broader homogenization (e.g., in an entire fuel block and reflector blocks with and without control rods) and the use of the usual approximations made to derive the theory. Therefore, it is prudent to assess the accuracy of the neutronics methods that are used for VHTR calculations. In this paper, we evaluate the accuracy of COMET for VHTR calculations using the prismatic VHTR benchmark problem (Connolly et al., 2015). COMET is a hybrid stochastic/deterministic transport method based on the incident flux response expansion theory (Mosher and Rahnema, 2006). The numerical implementation of the method (Zhang and Rahnema, 2012) consists of using response expansion coefficients (RC) generated by local calculations followed by core calculations to obtain the desired solution. A stochastic method is used to pre-compute a library of the incident flux RCs for all the unique coarse meshes (such as fuel and reflector blocks) in the core. Using this library, deterministic sweeps are performed to converge on the outgoing/incoming fluxes on coarse mesh interfaces in the core. The desired whole core solution is then generated by a superposition of the converged incident responses and the precomputed RCs. In the original method (Zhang and Rahnema, 2012), the phase space was treated continuously except for the energy variable. Extensive benchmarking of COMET
Corresponding author. E-mail address: [email protected] (F. Rahnema).
https://doi.org/10.1016/j.nucengdes.2019.110395 Received 11 July 2019; Received in revised form 15 October 2019; Accepted 21 October 2019 0029-5493/ © 2019 Published by Elsevier B.V.
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
(Zhang and Rahnema, 2017; Ulmer et al., 2016; Connolly and Rahnema, 2013) against the Monte Carlo method in various 3D problems typical of different reactor types has confirmed its high fidelity and formidable computational efficiency. Recently, COMET was extended (Rahnema and Zhang, 2018) to also treat the energy variable continuously. Like stochastic methods, this extension allows COMET to fully model the material heterogeneity down to the fuel and burnable absorber particle levels. As a result, in addition to any reactor type, it is demonstrated in this paper that COMET is an ideal tool for modeling the VHTRs because of its continuous-energy stochastic-like fidelity while retaining the formidable computational efficiency. The paper is organized as follows. For the sake of completeness and the reader’s convenience, the continuous energy COMET method is briefly reviewed in Section 2 and the VHTR benchmark problem is summarized in Section 3. The numerical results are presented and discussed in Section 4 and the paper is concluded in Section 5.
Rism ( r ,
Jis+, m =
1 F (r , k i
, E ) for r
Vi ,
,
, E) =
i
(r ,
Rismsm =
,
Jis±, m W (E )
, E) =
m( r
,
, E ) for r
m
m( r
(3)
In the above equation, represent the outgoing or incoming angular flux on surface s of coarse mesh i, denoted as Vis , m represents a set of orthogonal expansion functions, W (E ) is an asymptotic spectrum (Rahnema and Zhang, 2018), and Jis± , m are the corresponding expansion coefficients. Based on the incident flux response expansion theory, the flux distribution in each coarse mesh can be computed as the following superposition. ± is
i(r
,
Jis , m Rism ( r ,
, E) =
, E )for r
s, m
Vi .
, E) =
1 FRism ( r , k
, E)
s'm'
, m'
'
R m' m , is s
(7)
Vis
dr
nis+·
>0
d
dE (nis+· ) Rism ( r ,
, E)
m( r
,
, E)
,
, E ) = Pl (x ) Pn (y ) Up (µ ) Pq ( ) Bg,0 (E ),
(9)
3. Stylized prismatic VHTR benchmark problem A brief description of the stylized VHTR benchmark problem is presented in this section. Readers are referred to Connolly et al. (2015) for the detailed description. As shown in Fig. 1, the VHTR benchmark problem consists of a graphite inner reflector, three-rings of fuel blocks (yellow regions in Fig. 1), and an outer graphite reflector. The green regions in the reflectors represent the locations of the control channels, while the corresponding regions in the fuel blocks denote the reserved shutdown channels. In the axial direction, the core consists of a two-layer bottom reflector, a ten-layer active core, and a two-layer bottom reflector. The height of each layer is 79.3 cm. Vacuum boundary condition is imposed on all external surfaces. The following two core configurations are modeled in this paper. We note that a critical core configuration was not provided for the continuous energy case.
(4)
Here, Rism are the surface-to-volume response expansion functions that can be obtained by solving the following local fixed-source problem:
HRism ( r ,
is'
where the functions P and U are the Legendre polynomials and the Chebyshev polynomials of the second kind, respectively, and B represents the B-splines. That is, Legendre polynomials are used to expand the spatial dependence and a tensor product of Chebyshev polynomials of the second kind and Legendre polynomials are used to expand the polar and azimuthal angular distribution, respectively. It can be seen from Eqs. (3) and (9) that the energy distribution is expanded in terms of the products of the asymptotic spectrum W(E) and the 0th order B-splines with coarse energy bins. The numerical implementation of the COMET method consists of the following steps. In the first step, a continuous energy stochastic method is used to solve the local fixed-source problems defined by Eqs. (5) and (6) to obtain the library of the incident response coefficients. Because of the geometric robustness and the continuous energy treatment even the heterogeneity within the fuel and burnable absorber particles can be fully modeled. This detailed treatment of material heterogeneity is performed through the stochastic precomputation of the response coefficient (function) library. The precomputation of the library is similar to the industry practice of precomputation of the multigroup cross section library. In the second step, an iterative deterministic sweeping solver (Zhang and Rahnema, 2012) is used to converge the core solution via Eq. (7). In the last step, Eq. (4) is used to construct the flux (or the pin power) distribution in each coarse mesh by superposition.
(2)
Vis
J
In theory, any complete function set can be used in Eq. (3) to expand the outgoing/incoming angular fluxes on coarse mesh surfaces as long as they satisfy the orthogonal conditions. However, extensive benchmarks (Zhang and Rahnema, 2017; Ulmer et al., 2016; Connolly and Rahnema, 2013) have shown the following expansion functions are a good choice since they can sufficiently represent the angular flux phase distributions at a relatively low order and consequently achieve high computational efficiency while maintaining high accuracy.
where H and F are the standard transport and fission operators, k stands for the core eigenvalue (keff), i represents the angular flux in coarse mesh i at location r or (x , y , z ) in the direction of the unit vector or (µ, ) and with energy E , i is the flux from the neighboring meshes incident on the boundary of mesh i ( Vi ), and ni+ denotes the outward normal to Vi . The variables µ and are the cosine of the polar angle and the azimuthal angle, respectively. Since the incident flux i on each surface is not known a priori, it is assumed it can be expanded in terms of a set of orthogonal expansion functions as below. ± is ( r
(6)
Vis & · nis+ < 0
(8)
(1)
Vi & ·ni+ < 0,
, E )for r
, E )for r
where R m' m is the surface-to-surface response coefficients that can be is s computed by projecting the sourface-to-volume response functions m Ris ( r , , E ) onto the expansion basis.
with the boundary condition i(r
,
'
The detailed description of the COMET method and the continuous energy treatment can be found in Mosher and Rahnema (2006), Zhang and Rahnema (2012) and Rahnema and Zhang (2018). A short description is presented in this section for the sake of completeness and the readers convenience. To solve the Boltzmann transport equation in a reactor core, the spatial domain of the problem, V , is divided into a number of nonoverlapping coarse meshes {Vi i = 1, 2, , N } , and the phase space distribution of the flux in each coarse mesh can be obtained by solving the following set of local problems.
, E) =
m( r
The partial current expansion moments Jis+, m can be iteratively determined by projecting Eq. (2) onto the expansion basis as below.
2. Continuous energy COMET method
H i(r ,
, E ) = W (E )
(5)
with the boundary condition 2
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Fig. 1. Radial configuration of the VHTR benchmark problem; Taken from Connolly et al. (2015).
Fig. 3. Cut through, fuel pin; Taken from Connolly et al. (2015).
• All rods out (ARO): all control rods are withdrawn. • All rods in (ARI): all control rods in the inner and outer reflectors are
(assemblies) consist of small spheres (with a radius of 0.02125 cm) of uranium oxycarbide enriched to 15% in 235U by weight surrounded by four layers of materials: carbon, pyrolitic carbon, silicon carbide, and pyrolitic carbon. The packing fraction within the fuel pins is 35%. These fuel particles are modeled as a body centered cubic (BCC) lattice with the side length of 0.09664222 cm. it is noted that the choice of the BCC leads to unrealistic case of cutting particles/spheres as seen in Fig. 3. This of course is inconsequential from a benchmarking point of view.
fully inserted.
3.1. Fuel blocks The geometric configuration of the fuel blocks (with or without the reserve shutdown channel) is illustrated in Fig. 2. The apothem of each fuel block is 18 cm. Each fuel block consists of a hexagonal graphite structure (grey area in Fig. 2), 210 fuel pins (yellow regions), 108 coolant channels (white regions), 6 burnable absorber pins (green regions). Fuel blocks with the reserve shutdown channel have a similar geometric configuration as the standard fuel blocks except they have a 9.525-cm diameter hole located 9.756-cm away from the center. As a result, they have 186 fuel pins, 95 coolant channels, and 6 burnable absorber pins.
4. Numerical results The benchmark problem was modeled by both COMET and MCNP (X-5 Monte Carlo, 2005). The eigenvalues and the fission density distributions predicted by the two codes were then compared. The following four statistical quantities as defined in Zhang and Rahnema (2012) were used to summarize the differences between the pin fission density values (distributions) predicted by the two codes: the average (AVG), root mean square (RMS), mean relative/flux weighted (MRE), maximum (MAX) of the relative differences in the pin fission densities.
3.2. Fuel pins Fuel pins (see Fig. 3) imbedded in the graphite fuel blocks
Fig. 2. Configuration of fuel blocks and fuel blocks with the reserve shutdown channel; Taken from Connolly et al. (2015).
3
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Fig. 4. MCNP computation model.
4.1. MCNP calculation Because the core configurations have 1/12th symmetry in the radial direction and 1/2 symmetry in the axial direction, only 1/24th of the core was modeled in MCNP as shown in Fig. 4. This significantly reduced the MCNP computation time to maintain the same uncertainty as in the full core case. More importantly, this also avoided the MCNP fission source convergence problem that typically is a very challenging issue for large eigenvalue problems as fully discussed/analyzed in Zhang and Rahnema (2017). MCNP calculations were performed on an 84-CPU cluster. The fission source distribution was first converged by following 50 million particle histories, and then 25 million inactive and 500 million (4000 active cycles with 125,000 histories in each cycle) active particle histories were used to perform the actual calculations for each core configuration.
Fig. 5. Color set consisting of fuel and graphite blocks for generating asymptotic spectra.
surfaces. In the second approach, the spectrum on surface 5 from the color set calculation was used as the asymptotic function W (E ) on all fuel-fuel surfaces, the average of the incoming/outgoing spectra on surfaces 3 and 4 were used as the asymptotic functions on all fuel-reflector/reflector-fuel surfaces, and the average of the spectra on surfaces 1 and 2 was used as the asymptotic function on all reflector-reflector surfaces. The stochastic continuous energy response function generator (modified MCNP) was then used to pre-compute the response coefficients for all unique coarse meshes. The expansion order was set to {4,4,2,2} which represents the 4th order expansion on the two spatial variables on coarse mesh boundaries (surfaces) and the 2nd order expansion on the two angular (polar and azimuthal) variables. The energy distribution deviation away from W(E) was expanded by the 0th order B-splines with eight coarse energy bins with bin boundaries of 1.4572E−7, 6.2506E−7, 3.9279E−6, 1.3007E−4, 9.1188E−3, 8.2085E−1, 2.2313, and 20 MeV. It should be noted that the energy bin structure was originally used for the Pressurized Water Reactor benchmark problems (Rahnema et al., 2013) and may not be optimal for the VHTR problems. Regardless, it is demonstrated that the method works very well despite a choice that is not necessarily optimized or ideal. Of course, optimization of the energy bin structure is expected to further improve the COMET accuracy. However, this study is beyond the scope of this work. The response function library was precomputed with 15 million particle histories in each calculation that took about 3.2 h for a fuel block and 6 min for the other coarse meshes on a single CPU. The total computation time for the response function generation was about 873 h, or equivalently 10.4 h on the 84-CPU cluster.
4.2. COMET calculations Nine unique coarse meshes were used by COMET to model the VHTR benchmark problem. They include: two fuel blocks (without and with the reserve shutdown channel), two cooled reflector blocks (without and with the control rod hole), 3 hexagonal graphite reflector blocks (without the control rod hole, with an empty control rod hole, and with the control rod inserted), and two incomplete graphite blocks to match the core external boundary as shown in Fig. 1. The dimension and properties of each coarse mesh are given in Table 1. It can be seen from Section 2 that the choice of the asymptotic spectrum W (E ) is not unique. However, if W (E ) is very close to the actual spectrum, it can significantly lower the energy expansion order. In this work, the following two choices were tested: 1) a single asymptotic spectrum for all surfaces, and 2) surface-dependent asymptotic spectra based on its location in the core. A color set shown in Fig. 5 consisting of three graphite reflector blocks (coarse mesh 5 in Table 1) and three fuel blocks (coarse mesh 1 in Table 1) with the specular boundary condition was first used to generate the asymptotic spectra. In the first approach, the energy distribution on surface 5 (Fig. 5) was used as the asymptotic spectrum on all coarse mesh Table 1 Dimension and property of each unique coarse mesh. Index
Coarse mesh
Dimension: Apothem × Height (cm)
1 2 3 4 5 6 7 8 9
Fuel block Fuel block with reserve shutdown channel Cooled reflector Cooled reflector with reserve shutdown channel Reflector block Reflector block with control rod hole Reflector block with control rod inserted Side reflector Corner reflector
18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 18 × 49.3 1/2 of 18 × 49.3 hexagon 1/3 of 18 × 49.3 hexagon
4
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Table 2 Comparison of the core eigenvalues and computation times. Core
ARO
ARI
Method
COMET**
MCNP*
Eigenvalue Difference (pcm) Run Time (hours)
1.06617 ± 0.00004
Eigenvalue Difference (pcm) Run Time (hours)
0.91534 ± 0.00004
49.2
49.5
Single Asymptotic Spectrum
Surface-dependent Asymptotic Spectrum
1.06464 ± 0.00005 153 ± 6 0.89
1.06683 ± 0.00005 66 ± 6 0.90
0.91674 ± 0.00005 140 ± 6 0.85
0.91603 ± 0.00005 69 ± 6 0.85
* MCNP calculations modeled only 1/24th of the core on an 84-CPU cluster. It is estimated that a single full core calculation on this cluster would take about 50 days. ** COMET calculations modeled the full core on a single CPU of the cluster.
The response function library was used to perform core calculations for the two core configurations. To demonstrate the COMET computational efficiency, the full core was modeled in COMET. The core eigenvalue and the fuel pin fission density distributions computed by COMET are compared to the corresponding MCNP results in Section 4.3.
plotted in Figs. 6 and 7. Because of the axial symmetry only the top five axial layers are presented. For the reader’s interest, the pin fission density distributions predicted by COMET for the ARO and ARI cases are plotted in the appendix (Figs. A1 and A2), respectively. It can be seen from Figs. 6 and 7 that the maximum differences in the pin fission densities occurs in the peripheral fuel pins in the top layer next to the inner and outer reflectors when the single asymptotic spectrum was used. This is because a single asymptotic spectrum W (E ) with coarse energy bins is not sufficient for representing all surface spectra such as those on the fuel-fuel, fuel-reflector, reflector-reflector, reflector-vacuum, etc. Such an approximation leads to an underestimation of the thermal neutron flux on interfaces between the graphite blocks and the fuel blocks. As seen from Figs. 6b and 7b, the agreement between COMET and MCNP is excellent when the surface dependent asymptotic spectra were used on each surface without compromising the formidable computational efficiency of COMET. Note that the response function pre-computation time increased by a factor 2.
4.3. Comparison of calculated results The core eigenvalues predicted by COMET and MCNP and the corresponding computation times are presented in Table 2. It can be seen that the differences between the eigenvalues predicted by COMET with the single asymptotic spectrum and MCNP are 153 pcm and 140 pcm for the ARO and ARI cases, respectively. The corresponding differences for the surface-dependent spectra are 66 pcm and 69 pcm. The COMET full core calculations took only about 0.9 h on a single core, while the MCNP calculation for both of the 1/24th-core configurations took about 49 h each on the 84-CPU cluster as seen in Table 2. The comparison of the pin fission density distribution predicted by COMET and MCNP are statistically summarized in Table 3. The four statistical parameters are defined as follows: the average (AVG), the mean relative (MRE), the root mean square (RMS) and the maximum (MAX) relative differences. The average and maximum relative differences between COMET with the single asymptotic spectrum and MCNP are 1.04% and 3.72% for the ARO case and 0.98% and 4.60% for the ARI case, respectively. When the surface-dependent asymptotic spectra are used then the agreement between COMET and MCNP is improved significantly as expected. The corresponding differences are 0.36% and 1.84% for the ARO case and 0.53% and 2.13% for the ARI case. The difference between the COMET and MCNP pin fission density distributions are
5. Summary and conclusions In this paper, the performance of the recently developed continuous energy COMET method was evaluated for modeling VHTR problems. This was done by comparing the COMET solution (core eigenvalue and detailed fuel pin fission density distribution) to the corresponding MCNP solution in the stylized prismatic VHTR benchmark problem with control rods fully inserted and withdrawn. COMET treats the energy variable continuously by decomposing the angular flux incident on the coarse mesh surfaces as a product of a fast-varying function in energy, called the asymptotic spectrum, and a slowly varying function in energy. It was found that the use of a (color) set of asymptotic spectra (namely, fuel-fuel, fuel-reflector, reflector-fuel, and reflector-reflector) leads to excellent results. The differences between the eigenvalues calculated by COMET and MCNP were 66 pcm for the ARO case and 69 pcm for the ARI case. The average and maximum relative difference in the pin fission density distributions were 0.36% and 1.84% for the ARO case and 0.53% and 2.13% for the ARI case, respectively. COMET was also found to be significantly faster than MCNP. Each COMET full core calculation took about 0.9 h on a single CPU, while MCNP calculations for 1/24th of the core on an 84-CPU cluster took about 49 h each. It should be noted that if the full-core problems were to be modeled then at least 24 times more particle histories would have been needed to maintain the same statistical uncertainty. However, obtaining a fully symmetric (converged) solution will require averaging solutions from multiple (50–100) independent runs (Zhang and Rahnema, 2017). Conservatively, assuming a scaling factor of 1, this leads to about 6.8 years – 13.6 years of MCNP computation time on the 84-CPU cluster. This indicates that COMET is 4–5 orders of magnitude faster than MCNP. It can be concluded that COMET is an ideal tool for
Table 3 The statistics of the differences between the Pin fission density distribution predicted by COMET* and MCNP**. Core
ARO
ARI
Asymptotic Spectrum
Single
Surface-dependent
Single
Surface-dependent
AVG MRE RMS MAX
1.04% 1.09% 1.23% 3.72%
0.36% 0.34% 0.46% 1.84%
0.98% 0.95% 1.20% 4.60%
0.53% 0.52% 0.64% 2.13%
* The COMET pin fission density uncertainties are in the range of 0.042% – 0.064%. ** The MCNP pin fission density uncertainties are in the range of 0.16% – 0.30%. 5
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Fig. 6. Relative difference between the pin fission density distributions in the top five axial layers predicted by COMET and MCNP for the ARO case.
Fig. 7. Relative difference between the pin fission density distribution in the top five axial layers predicted by COMET and MCNP for the ARI case.
modeling the VHTRs because of its continuous-energy stochastic-like fidelity but with a formidable computational speed.
Acknowledgements The second author owns equity in a company that has licensed the COMET technologies from Georgia Tech. This study which is a demonstration of COMET could affect his personal financial status. The terms of this arrangement have been reviewed and approved by Georgia Tech in accordance with its conflict of interest policies.
Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 6
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
Appendix A. . Comet pin fission density distributions The normalized pin fission density distributions in the top five axial layers predicted by COMET for the ARO and ARI cases are plotted in Figs. A1 and A2, respectively.
Fig. A1. Normalized pin fission density distribution in the top five layers predicted by COMET with surface-dependent asymptotic spectra for the ARO case.
Fig. A2. Normalized pin fission density distribution in the top five layers predicted by COMET with surface-dependent asymptotic spectra for the ARI case.
7
Nuclear Engineering and Design 358 (2020) 110395
D. Zhang and F. Rahnema
References
reactor benchmark problem with UO2 and MOX fuel. Nucl. Technol. 184, 1–28. Ulmer, R., Rahnema, F., Connolly, K., 2016. A neutronic benchmark specification and COMET solution for the Advanced Burner Test Reactor. Ann. Nucl. Energy 87, 76–106. X-5 Monte Carlo Team, 2005. MCNP—A General Monte Carlo N-Particle Transport Code, Version 5. Los Alamos National Laboratory. Zhang, D., Rahnema, F., 2012. An efficient hybrid stochastic/deterministic coarse mesh neutron transport method. Ann. Nucl. Energy 41, 1–11. Zhang, D., Rahnema, F., 2017. Monte Carlo solution convergence analysis via COMET in a stylized integral inherently safe LWR benchmark problem. Ann. Nucl. Energy 100, 3–11.
Connolly, K., Rahnema, F., 2013. A heterogeneous coarse mesh radiation transport method for neutronic analysis of prismatic reactors. Ann. Nucl. Energy 56, 87–101. Connolly, K., Rahnema, F., Tsvetkov, P., 2015. Prismatic VHTR neutronic benchmark problems. Nucl. Eng. Des. 285, 207–240. Mosher, S., Rahnema, F., 2006. The incident flux response expansion method for heterogeneous coarse mesh transport problems. Transp. Theory Stat. Phys. 34, 1–26. Rahnema, F., Zhang, D., 2018. Continuous energy coarse mesh transport (COMET) method. Ann. Nuclear Energy 115. Rahnema, F., Hon, R., Douglass, S., 2013. A stylized three-dimensional pressurized water
8