- Email: [email protected]
BEAVRS full core burnup calculation in hot full power condition by RMC code
Annals of Nuclear Energy 101 (2017) 434–446
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
BEAVRS full core burnup calculation in hot full power condition by RMC code Shichang Liu a, Jingang Liang b, Qu Wu a, JuanJuan Guo a, Shanfang Huang a,⇑, Xiao Tang a, Zeguang Li c, Kan Wang a a b c
Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Nuclear Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China
a r t i c l e
i n f o
Article history: Received 2 July 2016 Received in revised form 20 November 2016 Accepted 21 November 2016 Available online 9 December 2016 Keywords: BEAVRS Neutronics/thermal-hydraulics coupling Sub-channel Large-scale detailed burnup calculation On-the-fly cross section treatment RMC
a b s t r a c t Monte Carlo method can provide high fidelity neutronics analysis of different types of nuclear reactors, owing to its advantages of the flexible geometry modeling and the use of continuous-energy nuclear cross sections. However, nuclear reactors are complex systems with multi-physics interacting and coupling. MC codes can couple with depletion solver and thermal-hydraulics (T/H) codes simultaneously for the ‘‘transport-burnup-thermal-hydraulics” coupling calculations. MIT BEAVRS is a typical ‘‘trans port-burnup-thermal-hydraulics” coupling benchmark. In this paper, RMC was coupled with subchannel code COBRA, equipped with on-the-fly temperature-dependent cross section treatment and large-scale detailed burnup calculation based on domain decomposition. Then RMC was applied to the full core burnup calculations of BEAVRS benchmark in hot full power (HFP) condition. The numerical tests show that domain decomposition method can achieve the consistent results compared with original version of RMC while enlarging the computational burnup regions. The results of HFP by RMC agree well with the reference values of BEAVRS benchmark and also agree well with those of MC21. This work proves the feasibility and accuracy of RMC in multi-physics coupling and lifecycle simulations of nuclear reactors. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The Monte Carlo method, which is often taken as a benchmark method to validate deterministic methods, is receiving rising attention due to its irreplaceable advantages such as flexibility in geometry treatment, the ability to use continuous-energy pointwise cross-sections, the easiness to parallelize and high-fidelity of simulations. Nuclear reactors are complex systems with multi-physics interacting and coupling. For examples, nuclides are generated or depleted during the lifecycle of reactors, and thermal-hydraulics has feedbacks on material temperature and density and thus nuclear cross sections. MC codes must be coupled with different solvers for high fidelity multi-physics simulations. In order to perform burnup calculations, MC codes are externally or internally coupled with a point burnup solver to figure out the nuclide densities in fuel during the depletion. It can also be coupled with ⇑ Corresponding author at: LIUQING Building Room 901, Tsinghua University, Beijing 100084, China. E-mail address: [email protected] (S. Huang). http://dx.doi.org/10.1016/j.anucene.2016.11.033 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.
thermal-hydraulics codes to obtain thermal-hydraulic feedback (Liu et al., 2015). When the burnup and TH feedbacks were taken into account simultaneously, it is known as ‘‘transport-burnup-th ermal-hydraulics” multi-physics coupling, which is crucial for realistic reactors simulations and benchmarks calculations such as MIT BEAVRS benchmark (Horelik and Herman, 2012). To consider the temperature effects on nuclear cross sections in high fidelity coupling calculations, recently, on-the-fly (OTF) techniques have been proposed in order to reduce the memory usage for both resolved resonance energy (Cullen and Weisbin, 1976; Yesilyurt et al., 2012; Yang et al., 2015; Forget et al., 2014) and thermal energy (Pavlou and Ji, 2014). To consider the distributed temperatures and coolant densities in the whole core, external and internal couplings are two traditional methods for neutronics/thermal-hydraulics (N-TH) coupling. External coupling is easily achieved but versatility is bad, however, internal coupling has good versatility but is more complex and needs a lot of code changes. Moreover, for three-dimensional pin-wise full core burnup calculations, the excessive memory demand is becoming a key obstacle for MC codes. To handle the memory problem, data parallel methods which include tally data decomposition (TDD)
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
(Liang et al., 2014) and spatial domain decomposition (SDD) (Liang et al., 2016) have been proposed. In this paper, a new hybrid coupling method has been developed to couple continuous-energy Monte Carlo code RMC (Wang et al., 2015) and the sub-channel code COBRA. In order to deal with the temperature dependence of cross sections in RMC code, the onthe-fly cross sections treatment has been developed for cross sections in resolved resonance region and thermal energy region. The spatial domain decomposition method based on nested MPI parallelism was developed to handle the memory problem of large-scale detailed full core burnup calculation. The ‘‘transport-burnup-ther mal-hydraulics” coupling system was applied to burnup simulations of BEAVRS benchmark in hot full power condition. The numerical tests show that domain decomposition can achieve the consistent results compared with original version. The results of HFP agree well with the reference values of BEAVRS benchmark and also agree well with those of MC21. This work proves the feasibility and accuracy of RMC code in multi-physics coupling and lifecycle simulations of nuclear reactors. The remainder of this paper is organized as follows. Section 2 introduces the methodology, including on-the-fly cross sections treatment for resolved resonance region (RRR) and thermal energy region, coupling codes, the proposed hybrid coupling method, capability of varying materials during burnup calculations and spatial domain decomposition based on nested MPI parallelism. In Section 3, the accuracy and efficacy of spatial domain decomposition was tested with burnup calculations of a PWR assembly case and a two-dimensional PWR full core case. Then the modeling of BEAVRS benchmark was established and the coupling system was applied to BEAVRS full core burnup calculations in hot full power condition. Finally, the conclusions are presented in Section 4.
2. Computational methods 2.1. On-the-fly temperature-dependent cross section treatment For thermal reactors such as PWR and HTGR, the temperature effects on nuclear cross sections at resolved resonance region and thermal energy region are significant. For resolved resonance region, Target motion sampling (TMS) method is a new on-thefly temperature treatment which takes the thermal motion of target nuclei into account explicitly to model arbitrary temperatures with only 0 K continuous-energy cross sections. TMS method (Viitanen and Leppänen, 2012) was firstly developed in Serpent code based on Woodcock tracking which is the main tracking routine in Serpent. TMS method was modified and developed in RMC code based on the ray tracking (Liu et al., 2016a). Moreover, TMS method was applied to reaction rate tally and depletion calculation for power generation feedback and nuclides densities evolutions. With TMS method, the reaction rates can be tallied with temperature feedback (Liu et al., 2016b). Temperature feedback in each burnup step can therefore be considered, so as to make the high fidelity ‘ ‘transport-burnup-thermal-hydraulics” coupling calculation. At thermal energies, the on-the-fly interpolation of thermal scattering data was developed in RMC to consider the thermal scattering and bound effect (Liu et al., 2016b). The S(a,b) data at the whole temperatures ranges are no longer needed to store, except for S(a,b) data at some reference temperatures. With TMS method and on-the-fly interpolation of thermal scattering, RMC has the capability of high fidelity ‘‘transport-burnupthermal-hydraulics” coupling for thermal reactors.
435
2.2. Coupling codes 2.2.1. Monte Carlo code RMC RMC is a continuous-energy Reactor Monte Carlo neutron and photon transport code being developed by Department of Engineering Physics at Tsinghua University, Beijing. The code RMC intends to solve reactor analysis problems, and is able to deal with complex geometry, using continuous energy point-wise cross sections of different materials and temperatures. As one of new generation Monte Carlo codes, RMC is aimed at achieving full core calculations and analysis with high fidelity and efficiency by means of advanced methodologies and algorithms as well as high performance computing techniques. RMC uses ACE format nuclear cross sections with S(a,b) and probability tables treatment in thermal and unresolved resonance energy ranges, adopts constructive solid geometry technique for flexible geometry modeling and employs ray-tracking method as main option for particle transport. 2.2.2. Thermal-hydraulics code COBRA COBRA is a sub-channel analysis code which computes the flow and enthalpy distributions in nuclear fuel rod bundles or cores for both steady state and transient conditions. COBRA has two subchannel models, for the PWR fuel assembly and PWR core analysis. RMC has been coupled with COBRA with fuel assembly model (Liu et al., 2015). In this paper, the PWR core model was used in which a whole fuel assembly and the associated coolant flow are treated as a single lumped-parameter channel. 2.3. Hybrid coupling method Consider the advantages and disadvantages of external and internal couplings, a new hybrid coupling method (Guo et al., 2016) is developed. Hybrid coupling means transforming data via external files of thermal hydraulics code and managing all the useful data by internal memory in neutronics code, as is shown in Fig. 1. RMC is based on ray tracking, in which neutron will stop at the boundaries between different cells. In N-TH coupling, different cell have different temperatures, and those cell filled by coolant have different densities. The fission power are generated in fuel pins which are surrounded by coolant, making up each individual pin cell. The pin cell will be further divided into axial meshes. Therefore, the whole geometry can be divided into three-dimensional meshes, and each mesh is an axial segment of each pin cell. The same mesh based geometry can be applied to both neutronics and thermal-hydraulics codes. In this way, four three-dimensional matrixes are stored in RMC which are for fuel temperatures, coolant temperatures, coolant densities and fission power, respectively. The values of the matrix are stored according to their coordinates in the whole geometry. Two flags of cells are used to distinguish fuel and coolant. The user-defined temperatures are used for cells with no flags.
Fig. 1. Schematic diagram of hybrid coupling.
436
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
During the transport routine, if neutrons are located in cells with fuel or coolant flags, the temperature and density of material it is located will be search in the temperature and density matrixes based on the indexes which are calculated according to the coordinates. Mesh tally is used in RMC to tally the fission power distributions. The coupling process is shown in Fig. 2 and can be described as follows: 1. The coupling program starts with initial power distribution; 2. The power distribution is written in the input file of COBRA, and the COBRA calculation is performed. Fuel temperature is computed in each rod, and coolant properties are calculated in the sub-channels; 3. When the COBRA run is finished, the output is read and the temperatures and densities are stored as three threedimensional matrixes in RMC; 4. Then RMC is run with three new matrixes of temperatures and densities. During the transport routine, the temperatures matrixes are used for one-the-fly temperature treatment of nuclides cross sections, and the densities matrixes are used for calculating the macro cross section for sampling the free flight distances. The power distributions in each fuel rod are tallied; 5. Compare the newly calculated power distribution with the power distribution of previous step to judge the convergence status. If the convergence criterion is reached, output the distributions of power, temperature and density. Otherwise, go to step 2 until it is converged.
2.4. Burnup capability of RMC RMC is developed with an embedded depletion module aimed at performing burnup calculations of large-scale problems with high efficiency. Several measures have been taken to strengthen the burnup capabilities of RMC:
Initial power distribution
COBRA calculate temperature & density distribution Transfer the power distribution to COBRA
Transfer the temperature and density to RMC
1. A new depletion module called DEPTH is developed and implemented. The DEPTH module is able to handle detailed depletion chains containing thousands of isotopes at an extremely fast speed with accuracy owing to its advanced matrixexponential solvers including the rational approximation methods and the Laguerre polynomial approximation method. Numerical results of several depletion cases have shown that the DEPTH module is more accurate and efficient than the formerly embedded depletion module (She et al., 2013). 2. The Energy-Bin method and the Cell-Mapping method are implemented to speed up the transport calculations with large numbers of nuclides and tally cells (She et al., 2013). 3. Auto expansion of burnable cells and materials in lattice structure (She et al., 2013). 4. The batch tally method and the parallelized depletion module have been utilized to better handle cases with massive amounts of burnup regions in parallel calculations (She et al., 2013). 5. Data decomposition technologies of reducing memory requirements for large-scale parallel MC burnup calculations (She et al., 2014). Besides the above capabilities, two capabilities have been developed for BEAVRS full core burnup calculation, including capacity of varying materials during burnup calculations and domain decomposition based on nested MPI parallelism. 2.5. Varying materials during burnup calculations In lifecycle simulations of PWR, the boron concentration in coolant keep adjusting and changing to maintain the critical condition of reactors. Therefore, the capacity of varying materials was developed in RMC in order to change the boron concentration in each burnup step. 2.6. Domain decomposition based on nested MPI parallelism Different from data decomposition method, spatial domain decomposition (SDD) divides spatial geometry into domains, which are simulated separately by parallel processors, and particles crossing domains are communicated for continuing tracking, as shown in Fig. 3. Fig. 4 shows the flow of SDD. In consideration of efficiency, an asynchronous particle communication algorithm is designed and implemented in RMC code (Liang et al., 2016). Furthermore, domain decomposition method was coupled with Monte Carlo burnup process under a strategy of utilizing consistent domain partition in both transport and depletion module. In other words, burnup regions can be decomposed automatically according to geometry decomposition in transport. Through the coupling, geometry-based data including data of tally, material and isotope densities, i.e. main memory-consuming data types are all decomposed.
Run RMC & tally power distribution
N Converged? Y Print power, temperature and density distribution Fig. 2. Coupling flow chart.
Fig. 3. Schematic diagram of particles transport with domain decomposition.
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
437
Fig. 4. Flow chart of domain decomposition method.
Consider the load balance and frequency of MPI communication, a whole geometry should not be divided into too many and too small domains. Therefore, nested MPI parallelism was used in
domain decomposition for large scale parallelism, as shown in Fig. 5. The particles parallelism can be performed independently with different processes in each domain, and domain replication
Fig. 5. Diagram of nested MPI parallelism.
438
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
was used in each domain. The parallel efficiency of particles parallelism was obviously better than parallelism of domain, as each particle can be treated independently. 3. Results and analysis 3.1. Verification of domain decomposition Before applying the domain decomposition strategy to the large scale burnup calculations, the accuracy and efficacy of domain decomposition was verified with burnup calculations of a PWR assembly case and a two-dimensional PWR full core case. 3.1.1. PWR assembly case Typical 17 17 PWR assembly with reflective boundary (Fig. 6) was calculated, which has 264 fuel pins. Transport and burnup coupled calculation was performed for 264 burnup regions are divided into four radial domains with two processes each domain. 2000 neutrons per cycle are used for totally 300 cycles with 50 inactive cycles, resulting in the standard deviation of 80 pcm on Kinf. The results with and without domain decomposition were compared in Fig. 7. It can be found that the maximum difference of kinf is about 250 pcm, and most of the differences are within three times of standard deviations (240 pcm). Moreover, the trend of difference fluctuates up or down, which should be the statistical fluctuations by Monte Carlo method.
Fig. 7. Kinf of PWR assembly case.
3.1.2. PWR core case A two-dimensional PWR core (Fig. 8) was calculated, which has 50,952 fuel pins. Transport and burnup coupled calculation was performed for 50,952 burnup regions are divided into 8 radial domains with 15 processes each domain. For the calculation without domain decomposition, the domain replication and particle parallel were used. The same number of processes (120 processes) were used for fair comparisons of calculation efficiency. 500,000 neutrons per cycle are used for totally 500 cycles with 200 inactive cycles, resulting in the standard deviation of 4.3 pcm on Keff. The results with and without domain decomposition were compared in Fig. 9. It can be found that the maximum difference of keff is about 9 pcm, and most of the differences are within three times
Fig. 8. Two-dimensional PWR core.
Fig. 6. Typical PWR assembly.
Fig. 9. Keff of PWR core case.
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
of standard deviations (12.9 pcm). Moreover, the trend of difference fluctuates up or down, which should be the statistical fluctuations by Monte Carlo method. Moreover, calculations times were compared in Fig. 10. The ratio of time between with and without domain decomposition increases with burnup, raising to 1.81. The results of these two cases show that domain decomposition can achieve consistent kinf/keff compared with original version, and the ratio of time between with and without domain decomposition increases with burnup. From Figs. 7 and 9, it can be found that the results with and without domain decomposition were not completely the same. The reason is that the random seed numbers of each particle will change due to domain decomposition. The same results could be attained by rearranging the random seed numbers of each particle. However, the rearrangement of random seed numbers has not been adopted in this paper. Although the results of with and without domain decomposition were not exactly the same, the results are within the range of statistical deviations of Monte Carlo methods. Therefore, the results with domain decomposition were consistent compared with original version of RMC.
439
Fig. 11. Radial cross section of BEAVRS core.
3.3. Results of HZP 3.2. BEAVRS benchmark modeling BEAVRS benchmark is specifications and measured results for two cycles of a PWR which was proposed by the MIT Computational Reactor Physics Group (Horelik and Herman, 2012). BEAVRS benchmark provides detailed descriptions and measured results for PWR with 193 fuel assemblies, and each 17 17 assembly contains 264 fuel rods, with three different enrichments. Both the measured data of hot zero power (HZP) and hot full power (HFP) are given in BEAVRS. The configuration of BEAVRS core is shown in Figs. 11and 12. And the parameters of core and assembly are listed in Tables 1 and 2. The detailed geometry of BEAVRS full core was built up by RMC, as shown in Figs. 13and 14. To consider the axial distributions of fission power and coolant density, the active core is divided into 10 axial segments. For thermal-hydraulics model in COBRA, only an octant of core has been considered for the 193 fuel assemblies in the full core. Therefore, the grid of 31 channels, each containing an individual fuel assembly, making up the lower triangular region of the quarter-core shown in Fig. 15. Each assembly is divided into 10 axial segments.
The HZP condition of BEARS was firstly calculated to verify the neutronics model and on-the-fly temperature treatment. The onthe-fly treatment was compared with the pre-generated cross section (XS) of exact temperature. For HZP, the coolant temperature is 566 K, the temperatures of other material (including fuel) are 293.15 K, and the boron concentration is 975 pcm. When using the pre-generated XS, XS of 293.15 K and 566 K should be prepared in advance, while only cross sections of 0 K are needed for on-thefly treatment. The results are compared in Table 3. It can be found that the Keff of on-the-fly method agrees well with the pregenerated XS. Keff of RMC also agrees well with the result of OpenMC. The relative radial pin power distributions were also compared with the results of MC21 (Kelly et al., 2013). The pin power distributions also have good agreement as shown in Fig. 16. For pre-generated XS, 200,000 neutrons per cycle are used for totally 400 cycles with 200 inactive cycles, while for on-the-fly method, 100,000 neutrons per cycle are used for totally 400 cycles with 200 inactive cycles. The ENDF/B7.0 cross-section library was used.
3.4. Results of HFP
Fig. 10. Calculations time of PWR core case.
The HFP condition of BEARS at zero burnup was calculated to verify the coupling between neutronics/thermal-hydraulics. For the coupling in HFP conditions, three feedback should be considered, including the temperature of coolant and fuel, the density of coolant and the boron concentration in coolant. For the boron concentration, there are two kinds of critical boron curves, one is in 100% core power and the other is with realistic power history, as is shown in Fig. 17. In the burnup calculations in HFP condition, the boron concentration was updated at every burnup step according to the given boron concentration of benchmark, to see whether the reactor is critical or not. However, for the HFP condition at zero burnup, there is no data of critical boron concentration, as the first burnup point in Fig. 17 for 100% power is for the fourth EFPD which is 599 pcm. Consider the effect of xenon poisoning at the beginning of burnup, the critical boron concentration must be larger than 599 pcm. Therefore, HFP condition of BEARS at zero burnup with different boron con-
440
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Fig. 12. Axial cross section BEAVRS core.
Table 1 Nominal value of thermal hydraulics condition. Item
Value
Units
Total power Initial mass flow rate Reference pressure Inlet water temperature Outlet water temperature
3411 17,083 155.17 292.78 310
MW kg/s bar °C °C
Table 2 Assembly size information. Item
Size
Number of fuel rods Number of guide tubes rods Active length (mm) Bundle pitch (mm) Fuel rod diameter (mm) Cladding inner diameter (mm) Cladding outer diameter (mm) Pin pitch (mm) Inner diameter of guide tube (mm) Outer diameter of guide tube (mm)
264 25 3657.6 215.04 7.84 8.00 9.14 12.60 11.22 12.04
Fig. 13. Radial cross section of BEAVRS core (RMC).
centrations were compared in Table 4. For HFP, 10,000 neutrons per cycle are used for totally 400 cycles with 150 inactive cycles. From Table 4, it can be found that with 975 pcm boron concentration, Keff of HZP is almost critical, while Keff of HFP is 1477.1 pcm smaller than 1. It shows the total effect of temperature feedback (including fuel, coolant and boron) is negative. On the other hand, 599 pcm is also not appropriate for HFP, as Keff is 3189.5 pcm larger. The reason is there is no xenon at zero burnup, while at the fourth EFPD the xenon poisoning reach saturation. By critical search of boron concentration, the critical boron concentration at zero burn is 850 pcm (see Fig. 18). In order to achieve the detailed power, 500,000 neutrons per cycle are used for totally 400 cycles with 150 inactive cycles. The convergence of different cycle neutrons number are compared in
Fig. 10. The ‘‘Max relative diff” (or ‘‘Max diff”) means the maximum relative difference of assembly power between two coupling iterations. It can be found that the ‘‘Max relative diff” of 500,000 neutrons is smaller. The statistic relative deviation of Monte Carlo calculations is also compared in Table 5. ‘‘Max RE” means the maximum relative deviation of assembly power in Monte Carlo calculations. It can be found that the convergence of coupling iterations depend on the statistic relative deviation of Monte Carlo calculations. When convergence is achieved, ‘‘Max diff” is within or close to three times of ‘‘Max RE”. The linear power per assembly, fuel average temperature, coolant temperature and coolant density at the middle of axial core are
441
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Fig. 17. Critical boron concentration evolution with burnup (reference). Fig. 14. Axial cross section of BEAVRS core (RMC). Table 4 Keff with different concentration.
Fig. 15. Radial channels of BEAVRS core (COBRA).
Table 3 Keff comparisons of pre-generated cross sections and on-the-fly method.
Pre-generated XS On-the-fly XS OpenMC
Keff
Diff
0.998207 (0.000103) 0.998455 (0.000178) 0.99920 (0.00004)
– 0.000248 0.000993
power
Condition
Keff
HZP + 975 ppm HFP + 975 ppm HFP + 599 ppm HFP + 850 ppm
0.998461 0.985229 1.031895 0.999847
and
boron
(0.000070) (0.000498) (0.000445) (0.000457)
shown in Figs. 19–22. It can be found that the fuel average temperature and coolant temperature have the checkerboard distributions similar to linear power distributions, proving that the results of thermal-hydraulics calculations are reasonable. The axial distributions of fuel average temperature, coolant temperature and coolant density in the hottest assembly are shown in Fig. 23. The detailed radial pin power distributions of RMC in HFP condition at zero burnup are compared with that of HZP and the results of MC21 (Daniel et al., 2014), as shown in Fig. 24. Compared with the HZP results, the power distributions of HFP are more uniform, due to the negative temperature feedback. The power distributions of RMC agree well with the results of MC21 although there are some differences. The differences are because the burnup of
RMC
MC21 Fig. 16. Relative radial pin power distributions (HZP).
442
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Another important phenomenon should be noted is that unnatural and unreasonable cross shaped power peaks appear in some hot assemblies in the middle of core. The reason is that in the thermal-hydraulics model of MC21, each assembly is divided into 4 separate channels, which gives raise to the power peaks between channels in the same assembly. 3.5. Results of burnup calculation in HFP
Fig. 18. Convergence curve of different cycle neutrons number.
Table 5 Convergence status of different cycle neutrons number. Cycle neutrons number
Max RE
Max diff
Diff/RE
500,000 10,000
0.0231 0.152
0.03139 0.26459
1.359 1.741
MC21 results is 92 calendar days (1023.3 MWD/THM) with 76% full power and 702.5 ppm boron concentration.
On the basic of neutronics/thermal-hydraulics coupling in HFP condition, the domain decomposition strategy was added to perform the full core ‘‘transport-burnup-thermal-hydraulics” of BEAVRS benchmark. Each fuel pin and burnable poison pin is treated as individual burnup regions, and each pin was divided into 10 axial segments, summing up to 640,000 burnup regions. The whole core is divided into 16 domains: 2 axial domains and 8 radial domains, as shown in Fig. 25. The large scale parallelism was performed on Tianhe2 super computer. 30 nodes with 720 processes were used, with 45 processes in each domain. For the power normalization of depletion calculation, the ratio of real power and relative power tallied by Monte Carlo code should be figured out. For domain decomposition, each domain has its total power, and communications between domains will be conducted with MPI_Allreduce to sum up the domain power to get the total tallied power. Once the total tallied power was attained, the ratio of real power and relative power can be figured out. It should be noted that in this MPI communications, only the domain power should be communicated, therefore the data size is small.
Fig. 19. Linear power per assembly (w/m).
Fig. 20. Fuel average temperature (K).
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
443
Fig. 21. Coolant temperature (K).
Fig. 22. Coolant density (g/cm3).
Fig. 23. Axial distributions.
The ‘‘transport-burnup-thermal-hydraulics” strategy was shown in Table 6. Two sequential transport/thermal-hydraulics coupling iterations were performed firstly, to attain the power distributions with thermal-hydraulics feedback considered, as shown in Fig. 26. Then the burnup calculation were performed, using the predictor-corrector depletion strategy. As the boron concentration in coolant changes in different burnup steps, the capacity of varying materials was used to change the boron concentration in each burnup step according to Fig. 17 to keep the reactor critical.
For the hot full power calculations, the specific power is 41.699 W/gHM for each burnup step. The burnup steps are: 4, 7, 5, 6, 9, 5, 16, 17, 16, 11, 14, 14 (days). . ., which is set according the reference of critical boron concentration in Fig. 16. This setting of burnup steps is convenient for the comparisons of reactivity with critical boron concentrations. To investigate the effects of thermal-hydraulics feedback in burnup calculations, another burnup calculation without thermalhydraulics feedback was performed, while changing the boron concentration according to Fig. 17. For the case without thermalhydraulics feedback, the fixed temperature cross sections were used, in which the coolant temperature is set as 600 K, and fuel temperature is 900 K. The results are compared in Fig. 27. It was found that keff of case without feedback keeps increasing when burnp increases, raising to 3% larger than critical. While keff with feedback always maintains within 0.99–1.0. It can be concluded that the critical boron concentration with thermal-hydraulics feedback by RMC agrees with the critical boron concentration given by the benchmark (Fig. 17). The difference of reactivity may be resulted from a lot of reasons. But the main reason is because of the calculations were carried on by hot full (100%) power. However, in the cycle 1 of BEAVRS, the real power is not constant but changed with burnup and the average power is about 75%. When increasing the power from 75% to 100%, keff will decrease because the temperature is higher and negative reactivity feedback. Moreover, when the power increase, the saturated xenon concentration will increase, also leading to the decrease of keff. The Relative radial pin power distributions with different burnup were calculated by RMC and compared in Fig. 28. It can be
444
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
HZP
HFP(RMC, day 0)
HFP(MC21, day 92)
Fig. 24. Relative radial pin power distributions (HFP).
Fig. 25. Diagram of domain decomposition of BEAVRS.
Table 6 Coupling strategy. Iteration
Contents
Calculation conditions
1
1st transport/1st thermalhydraulics 2nd transport /2nd thermalhydraulics Predictor depletion Corrector transport/corrector depletion
106 neutrons per cycle, 150 inactive, 100 active 106 neutrons per cycle, 150 inactive, 400 active
2 3 4
106 neutrons per cycle, 150 inactive, 400 active
found that the power distributions of core becomes uniform when the burnp increases. 4. Conclusions The‘‘transport-burnup-thermal-hydraulics” multi-physics coupling is crucial for realistic reactors simulations and benchmarks
calculations such as MIT BEAVRS benchmark. In order to perform the high fidelity‘‘transport-burnup-thermal-hydraulics” coupling, several advanced techniques were developed in RMC. The onthe-fly temperature-dependent cross section treatment was developed to consider the temperature effects on nuclear cross sections in resolved resonance region and thermal energy region. A hybrid coupling method was proposed and RMC was coupled with subchannel code COBRA to consider the distributed temperatures and coolant densities in the whole core. The hybrid coupling method can reduce the difficulty of modeling and improve the versatility of coupling by managing all the useful data by internal memory in neutronics code, while making good use of the existing thermal-hydraulics codes. The spatial domain decomposition method based on nested MPI parallelism was developed to handle the memory problem of pin-wise three-dimensional full core burnup calculations. The ‘‘transport-burnup-thermal-hydraulics” coupling system was applied to burnup simulations of BEAVRS benchmark in hot full power condition.
445
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Fig. 27. Keff evolution with burnup.
Fig. 26. Flowchart for the transport-T/H-burnup coupling scheme.
The proposed domain decomposition method based on nested MPI parallelism was firstly verified by solving the burnup problems of a typical PWR assembly and a two-dimensional PWR core. The results of these two cases show that domain decomposition can achieve consistent kinf/keff compared with original version of RMC, while enlarging the computational burnup regions. The time ratio of with and without domain decomposition increases with burnup. The on-the-fly temperature treatment was then verified by the hot zero power condition of BEAVRS. Keff of on-the-fly method agrees well with the pre-generated cross sections of exact temperature. The relative radial pin power distributions also agree with the results of MC21.
Burnup=0
For the HFP calculations, the detailed radial pin power distributions of RMC in HFP condition at zero burnup are compared with those of HZP and the results of MC21. The critical boron concentration of HFP at zero burnup was found to be 850pcm. The convergence of coupling iterations are compared with statistic relative deviation of Monte Carlo calculations, showing that the convergence of coupling iterations depends on the statistic relative deviation of Monte Carlo calculations. The power distributions of HFP are more uniform than those of HZP due to the negative temperature feedback. Moreover, unreasonable cross shaped power peaks are found in some hot assemblies of MC21 results, which are resulted from the thermal-hydraulics model of MC21. Finally, the burnup calculations in HFP were performed. The results show that keff with feedback always maintains within 0.99–1.0 when burnup increases, while keff of case without feedback keeps increasing when burnup increases, rising to 3% larger than criticality. It can be concluded that the critical boron concentration with thermalhydraulics feedback by RMC agrees with the critical boron concentration given by the benchmark. The relative radial pin power distributions with different burnup were also compared and found that the power distributions of core become uniform when the burnup increases. This work provides the foundation for ‘‘transport-burnup-ther mal-hydraulics” coupling calculations in multi-physics coupling and lifecycle simulations of nuclear reactors using RMC code. More
Burnup=1501 MWD/THM Fig. 28. Relative radial pin power distributions with different burnup.
Burnup=3544 MWD/THM
446
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
research will be done to fulfill the two cycle simulations of BEAVRS benchmark with RMC code, and more advanced thermal-hydraulic solver such as COBRA-TF would be used for coupling with pin-bypin sub-channel and parallelism in TH calculations. Acknowledgments This work is partially supported by National Science and Technology Major Project of Large Advanced PWR Nuclear Power Plant in China (2011ZX06004-024) and Project 91126017/11475098 by National Natural Science Foundation of China. References Cullen, D.E., Weisbin, C.R., 1976. Exact Doppler broadening of tabulated cross sections. Nucl. Sci. Eng. 60, 199. Daniel, K., Brian, A., Paul, R., 2014. Analysis of select BEAVRS PWR benchmark cycle 1 results using MC21 and OpenMC. In: PHYSOR2014, Kyoto, Japan, September 28–October 3, 2014. Forget, B., Xu, S., Smith, K., 2014. Direct Doppler broadening in Monte Carlo simulations using the multipole representation. Ann. Nucl. Energy 64, 78–85. Horelik, N., Herman, B., 2012. Benchmark for evaluation and validation of reactor simulations. MIT Computational Reactor Physics Group. URL:. Guo, J. et al., 2016. Neutronics/thermal-hydraulics coupling with RMC and CTF for BEAVRS benchmark calculation. In: 2016 ANS Winter Meeting. Las Vegas. Kelly, D.J. et al., 2013. MC21 analysis of the MIT PWR benchmark: hot zero power results. In: Proceedings of the International Conference on Mathematics and Com-Putational Methods Applied to Nuclear Science and Engineering. Sun Valley, Idaho. May 5–9.
Liang, J., Wang, K., Yu, G., et al., 2014. Research on tally data decomposition algorithms based on reactor Monte Carlo code RMC. Nucl. Power Eng. 35 (4), 142 (in Chinese). Liang, J., Wang, K., Qiu, Y., Chai, X., Qiang, S., 2016. Domain decomposition strategy for pin-wise full-core Monte Carlo depletion calculation with the reactor Monte Carlo code. Nucl. Eng. Technol. 48, 635–641. Liu, S.C., Yu, J.K., Liang, J.G., Wang, K., 2015. Study of neutronics and thermalhydraulics coupling with RMC COBRA-EN system. In: 7th International Conference on Modelling and Simulation in Nuclear Science and Engineering (7ICMSNSE), Ottawa, Ontario, Canada, October 18–21, 2015. Liu, S.C., Yuan, Y., Yu, J.K., Wang, K., 2016a. Development of on-the-fly temperaturedependent cross-sections treatment in RMC code. Ann. Nucl. Energy 94, 144– 149. Liu, S.C., Yuan, Y., Yu, J.K., Wang, K., 2016b. Reaction rate tally and depletion calculation with on-the-fly temperature treatment. Ann. Nucl. Energy 92, 277– 283. Pavlou, T.A., Ji, W., 2014. On-the-fly sampling of temperature-dependent thermal neutron scattering data for Monte Carlo simulations. Ann. Nucl. Energy 71, 411– 426. She, D. et al., 2013. Development of burnup methods and capabilities in Monte Carlo code RMC. Ann. Nucl. Energy 51, 289–294. She, D. et al., 2014. 2D full-core Monte Carlo pin-by-pin burnup calculations with the RMC code. Ann. Nucl. Energy 64, 201–205. Viitanen, T., Leppänen, J., 2012. Explicit treatment of thermal motion in continuousenergy Monte Carlo tracking routines. Nucl. Sci. Eng. 171, 165. Wang, K. et al., 2015. RMC – a Monte Carlo code for reactor core analysis. Ann. Nucl. Energy 82, 121–129. Yang, F., Liang, J.G., Yu, G.L., Wang, K., 2015. Research on on-the-fly Doppler broadening based on the stochastic algorithm. In: 2015 ANS Annual Meeting. San Antonio, TX. June 7–11. Yesilyurt, G., Martin, W.R., Brown, F., 2012. On-the-fly Doppler broadening for Monte Carlo codes. Nucl. Sci. Eng. 171, 239.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
BEAVRS full core burnup calculation in hot full power condition by RMC code Shichang Liu a, Jingang Liang b, Qu Wu a, JuanJuan Guo a, Shanfang Huang a,⇑, Xiao Tang a, Zeguang Li c, Kan Wang a a b c
Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Nuclear Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China
a r t i c l e
i n f o
Article history: Received 2 July 2016 Received in revised form 20 November 2016 Accepted 21 November 2016 Available online 9 December 2016 Keywords: BEAVRS Neutronics/thermal-hydraulics coupling Sub-channel Large-scale detailed burnup calculation On-the-fly cross section treatment RMC
a b s t r a c t Monte Carlo method can provide high fidelity neutronics analysis of different types of nuclear reactors, owing to its advantages of the flexible geometry modeling and the use of continuous-energy nuclear cross sections. However, nuclear reactors are complex systems with multi-physics interacting and coupling. MC codes can couple with depletion solver and thermal-hydraulics (T/H) codes simultaneously for the ‘‘transport-burnup-thermal-hydraulics” coupling calculations. MIT BEAVRS is a typical ‘‘trans port-burnup-thermal-hydraulics” coupling benchmark. In this paper, RMC was coupled with subchannel code COBRA, equipped with on-the-fly temperature-dependent cross section treatment and large-scale detailed burnup calculation based on domain decomposition. Then RMC was applied to the full core burnup calculations of BEAVRS benchmark in hot full power (HFP) condition. The numerical tests show that domain decomposition method can achieve the consistent results compared with original version of RMC while enlarging the computational burnup regions. The results of HFP by RMC agree well with the reference values of BEAVRS benchmark and also agree well with those of MC21. This work proves the feasibility and accuracy of RMC in multi-physics coupling and lifecycle simulations of nuclear reactors. Ó 2016 Elsevier Ltd. All rights reserved.
1. Introduction The Monte Carlo method, which is often taken as a benchmark method to validate deterministic methods, is receiving rising attention due to its irreplaceable advantages such as flexibility in geometry treatment, the ability to use continuous-energy pointwise cross-sections, the easiness to parallelize and high-fidelity of simulations. Nuclear reactors are complex systems with multi-physics interacting and coupling. For examples, nuclides are generated or depleted during the lifecycle of reactors, and thermal-hydraulics has feedbacks on material temperature and density and thus nuclear cross sections. MC codes must be coupled with different solvers for high fidelity multi-physics simulations. In order to perform burnup calculations, MC codes are externally or internally coupled with a point burnup solver to figure out the nuclide densities in fuel during the depletion. It can also be coupled with ⇑ Corresponding author at: LIUQING Building Room 901, Tsinghua University, Beijing 100084, China. E-mail address: [email protected] (S. Huang). http://dx.doi.org/10.1016/j.anucene.2016.11.033 0306-4549/Ó 2016 Elsevier Ltd. All rights reserved.
thermal-hydraulics codes to obtain thermal-hydraulic feedback (Liu et al., 2015). When the burnup and TH feedbacks were taken into account simultaneously, it is known as ‘‘transport-burnup-th ermal-hydraulics” multi-physics coupling, which is crucial for realistic reactors simulations and benchmarks calculations such as MIT BEAVRS benchmark (Horelik and Herman, 2012). To consider the temperature effects on nuclear cross sections in high fidelity coupling calculations, recently, on-the-fly (OTF) techniques have been proposed in order to reduce the memory usage for both resolved resonance energy (Cullen and Weisbin, 1976; Yesilyurt et al., 2012; Yang et al., 2015; Forget et al., 2014) and thermal energy (Pavlou and Ji, 2014). To consider the distributed temperatures and coolant densities in the whole core, external and internal couplings are two traditional methods for neutronics/thermal-hydraulics (N-TH) coupling. External coupling is easily achieved but versatility is bad, however, internal coupling has good versatility but is more complex and needs a lot of code changes. Moreover, for three-dimensional pin-wise full core burnup calculations, the excessive memory demand is becoming a key obstacle for MC codes. To handle the memory problem, data parallel methods which include tally data decomposition (TDD)
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
(Liang et al., 2014) and spatial domain decomposition (SDD) (Liang et al., 2016) have been proposed. In this paper, a new hybrid coupling method has been developed to couple continuous-energy Monte Carlo code RMC (Wang et al., 2015) and the sub-channel code COBRA. In order to deal with the temperature dependence of cross sections in RMC code, the onthe-fly cross sections treatment has been developed for cross sections in resolved resonance region and thermal energy region. The spatial domain decomposition method based on nested MPI parallelism was developed to handle the memory problem of large-scale detailed full core burnup calculation. The ‘‘transport-burnup-ther mal-hydraulics” coupling system was applied to burnup simulations of BEAVRS benchmark in hot full power condition. The numerical tests show that domain decomposition can achieve the consistent results compared with original version. The results of HFP agree well with the reference values of BEAVRS benchmark and also agree well with those of MC21. This work proves the feasibility and accuracy of RMC code in multi-physics coupling and lifecycle simulations of nuclear reactors. The remainder of this paper is organized as follows. Section 2 introduces the methodology, including on-the-fly cross sections treatment for resolved resonance region (RRR) and thermal energy region, coupling codes, the proposed hybrid coupling method, capability of varying materials during burnup calculations and spatial domain decomposition based on nested MPI parallelism. In Section 3, the accuracy and efficacy of spatial domain decomposition was tested with burnup calculations of a PWR assembly case and a two-dimensional PWR full core case. Then the modeling of BEAVRS benchmark was established and the coupling system was applied to BEAVRS full core burnup calculations in hot full power condition. Finally, the conclusions are presented in Section 4.
2. Computational methods 2.1. On-the-fly temperature-dependent cross section treatment For thermal reactors such as PWR and HTGR, the temperature effects on nuclear cross sections at resolved resonance region and thermal energy region are significant. For resolved resonance region, Target motion sampling (TMS) method is a new on-thefly temperature treatment which takes the thermal motion of target nuclei into account explicitly to model arbitrary temperatures with only 0 K continuous-energy cross sections. TMS method (Viitanen and Leppänen, 2012) was firstly developed in Serpent code based on Woodcock tracking which is the main tracking routine in Serpent. TMS method was modified and developed in RMC code based on the ray tracking (Liu et al., 2016a). Moreover, TMS method was applied to reaction rate tally and depletion calculation for power generation feedback and nuclides densities evolutions. With TMS method, the reaction rates can be tallied with temperature feedback (Liu et al., 2016b). Temperature feedback in each burnup step can therefore be considered, so as to make the high fidelity ‘ ‘transport-burnup-thermal-hydraulics” coupling calculation. At thermal energies, the on-the-fly interpolation of thermal scattering data was developed in RMC to consider the thermal scattering and bound effect (Liu et al., 2016b). The S(a,b) data at the whole temperatures ranges are no longer needed to store, except for S(a,b) data at some reference temperatures. With TMS method and on-the-fly interpolation of thermal scattering, RMC has the capability of high fidelity ‘‘transport-burnupthermal-hydraulics” coupling for thermal reactors.
435
2.2. Coupling codes 2.2.1. Monte Carlo code RMC RMC is a continuous-energy Reactor Monte Carlo neutron and photon transport code being developed by Department of Engineering Physics at Tsinghua University, Beijing. The code RMC intends to solve reactor analysis problems, and is able to deal with complex geometry, using continuous energy point-wise cross sections of different materials and temperatures. As one of new generation Monte Carlo codes, RMC is aimed at achieving full core calculations and analysis with high fidelity and efficiency by means of advanced methodologies and algorithms as well as high performance computing techniques. RMC uses ACE format nuclear cross sections with S(a,b) and probability tables treatment in thermal and unresolved resonance energy ranges, adopts constructive solid geometry technique for flexible geometry modeling and employs ray-tracking method as main option for particle transport. 2.2.2. Thermal-hydraulics code COBRA COBRA is a sub-channel analysis code which computes the flow and enthalpy distributions in nuclear fuel rod bundles or cores for both steady state and transient conditions. COBRA has two subchannel models, for the PWR fuel assembly and PWR core analysis. RMC has been coupled with COBRA with fuel assembly model (Liu et al., 2015). In this paper, the PWR core model was used in which a whole fuel assembly and the associated coolant flow are treated as a single lumped-parameter channel. 2.3. Hybrid coupling method Consider the advantages and disadvantages of external and internal couplings, a new hybrid coupling method (Guo et al., 2016) is developed. Hybrid coupling means transforming data via external files of thermal hydraulics code and managing all the useful data by internal memory in neutronics code, as is shown in Fig. 1. RMC is based on ray tracking, in which neutron will stop at the boundaries between different cells. In N-TH coupling, different cell have different temperatures, and those cell filled by coolant have different densities. The fission power are generated in fuel pins which are surrounded by coolant, making up each individual pin cell. The pin cell will be further divided into axial meshes. Therefore, the whole geometry can be divided into three-dimensional meshes, and each mesh is an axial segment of each pin cell. The same mesh based geometry can be applied to both neutronics and thermal-hydraulics codes. In this way, four three-dimensional matrixes are stored in RMC which are for fuel temperatures, coolant temperatures, coolant densities and fission power, respectively. The values of the matrix are stored according to their coordinates in the whole geometry. Two flags of cells are used to distinguish fuel and coolant. The user-defined temperatures are used for cells with no flags.
Fig. 1. Schematic diagram of hybrid coupling.
436
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
During the transport routine, if neutrons are located in cells with fuel or coolant flags, the temperature and density of material it is located will be search in the temperature and density matrixes based on the indexes which are calculated according to the coordinates. Mesh tally is used in RMC to tally the fission power distributions. The coupling process is shown in Fig. 2 and can be described as follows: 1. The coupling program starts with initial power distribution; 2. The power distribution is written in the input file of COBRA, and the COBRA calculation is performed. Fuel temperature is computed in each rod, and coolant properties are calculated in the sub-channels; 3. When the COBRA run is finished, the output is read and the temperatures and densities are stored as three threedimensional matrixes in RMC; 4. Then RMC is run with three new matrixes of temperatures and densities. During the transport routine, the temperatures matrixes are used for one-the-fly temperature treatment of nuclides cross sections, and the densities matrixes are used for calculating the macro cross section for sampling the free flight distances. The power distributions in each fuel rod are tallied; 5. Compare the newly calculated power distribution with the power distribution of previous step to judge the convergence status. If the convergence criterion is reached, output the distributions of power, temperature and density. Otherwise, go to step 2 until it is converged.
2.4. Burnup capability of RMC RMC is developed with an embedded depletion module aimed at performing burnup calculations of large-scale problems with high efficiency. Several measures have been taken to strengthen the burnup capabilities of RMC:
Initial power distribution
COBRA calculate temperature & density distribution Transfer the power distribution to COBRA
Transfer the temperature and density to RMC
1. A new depletion module called DEPTH is developed and implemented. The DEPTH module is able to handle detailed depletion chains containing thousands of isotopes at an extremely fast speed with accuracy owing to its advanced matrixexponential solvers including the rational approximation methods and the Laguerre polynomial approximation method. Numerical results of several depletion cases have shown that the DEPTH module is more accurate and efficient than the formerly embedded depletion module (She et al., 2013). 2. The Energy-Bin method and the Cell-Mapping method are implemented to speed up the transport calculations with large numbers of nuclides and tally cells (She et al., 2013). 3. Auto expansion of burnable cells and materials in lattice structure (She et al., 2013). 4. The batch tally method and the parallelized depletion module have been utilized to better handle cases with massive amounts of burnup regions in parallel calculations (She et al., 2013). 5. Data decomposition technologies of reducing memory requirements for large-scale parallel MC burnup calculations (She et al., 2014). Besides the above capabilities, two capabilities have been developed for BEAVRS full core burnup calculation, including capacity of varying materials during burnup calculations and domain decomposition based on nested MPI parallelism. 2.5. Varying materials during burnup calculations In lifecycle simulations of PWR, the boron concentration in coolant keep adjusting and changing to maintain the critical condition of reactors. Therefore, the capacity of varying materials was developed in RMC in order to change the boron concentration in each burnup step. 2.6. Domain decomposition based on nested MPI parallelism Different from data decomposition method, spatial domain decomposition (SDD) divides spatial geometry into domains, which are simulated separately by parallel processors, and particles crossing domains are communicated for continuing tracking, as shown in Fig. 3. Fig. 4 shows the flow of SDD. In consideration of efficiency, an asynchronous particle communication algorithm is designed and implemented in RMC code (Liang et al., 2016). Furthermore, domain decomposition method was coupled with Monte Carlo burnup process under a strategy of utilizing consistent domain partition in both transport and depletion module. In other words, burnup regions can be decomposed automatically according to geometry decomposition in transport. Through the coupling, geometry-based data including data of tally, material and isotope densities, i.e. main memory-consuming data types are all decomposed.
Run RMC & tally power distribution
N Converged? Y Print power, temperature and density distribution Fig. 2. Coupling flow chart.
Fig. 3. Schematic diagram of particles transport with domain decomposition.
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
437
Fig. 4. Flow chart of domain decomposition method.
Consider the load balance and frequency of MPI communication, a whole geometry should not be divided into too many and too small domains. Therefore, nested MPI parallelism was used in
domain decomposition for large scale parallelism, as shown in Fig. 5. The particles parallelism can be performed independently with different processes in each domain, and domain replication
Fig. 5. Diagram of nested MPI parallelism.
438
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
was used in each domain. The parallel efficiency of particles parallelism was obviously better than parallelism of domain, as each particle can be treated independently. 3. Results and analysis 3.1. Verification of domain decomposition Before applying the domain decomposition strategy to the large scale burnup calculations, the accuracy and efficacy of domain decomposition was verified with burnup calculations of a PWR assembly case and a two-dimensional PWR full core case. 3.1.1. PWR assembly case Typical 17 17 PWR assembly with reflective boundary (Fig. 6) was calculated, which has 264 fuel pins. Transport and burnup coupled calculation was performed for 264 burnup regions are divided into four radial domains with two processes each domain. 2000 neutrons per cycle are used for totally 300 cycles with 50 inactive cycles, resulting in the standard deviation of 80 pcm on Kinf. The results with and without domain decomposition were compared in Fig. 7. It can be found that the maximum difference of kinf is about 250 pcm, and most of the differences are within three times of standard deviations (240 pcm). Moreover, the trend of difference fluctuates up or down, which should be the statistical fluctuations by Monte Carlo method.
Fig. 7. Kinf of PWR assembly case.
3.1.2. PWR core case A two-dimensional PWR core (Fig. 8) was calculated, which has 50,952 fuel pins. Transport and burnup coupled calculation was performed for 50,952 burnup regions are divided into 8 radial domains with 15 processes each domain. For the calculation without domain decomposition, the domain replication and particle parallel were used. The same number of processes (120 processes) were used for fair comparisons of calculation efficiency. 500,000 neutrons per cycle are used for totally 500 cycles with 200 inactive cycles, resulting in the standard deviation of 4.3 pcm on Keff. The results with and without domain decomposition were compared in Fig. 9. It can be found that the maximum difference of keff is about 9 pcm, and most of the differences are within three times
Fig. 8. Two-dimensional PWR core.
Fig. 6. Typical PWR assembly.
Fig. 9. Keff of PWR core case.
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
of standard deviations (12.9 pcm). Moreover, the trend of difference fluctuates up or down, which should be the statistical fluctuations by Monte Carlo method. Moreover, calculations times were compared in Fig. 10. The ratio of time between with and without domain decomposition increases with burnup, raising to 1.81. The results of these two cases show that domain decomposition can achieve consistent kinf/keff compared with original version, and the ratio of time between with and without domain decomposition increases with burnup. From Figs. 7 and 9, it can be found that the results with and without domain decomposition were not completely the same. The reason is that the random seed numbers of each particle will change due to domain decomposition. The same results could be attained by rearranging the random seed numbers of each particle. However, the rearrangement of random seed numbers has not been adopted in this paper. Although the results of with and without domain decomposition were not exactly the same, the results are within the range of statistical deviations of Monte Carlo methods. Therefore, the results with domain decomposition were consistent compared with original version of RMC.
439
Fig. 11. Radial cross section of BEAVRS core.
3.3. Results of HZP 3.2. BEAVRS benchmark modeling BEAVRS benchmark is specifications and measured results for two cycles of a PWR which was proposed by the MIT Computational Reactor Physics Group (Horelik and Herman, 2012). BEAVRS benchmark provides detailed descriptions and measured results for PWR with 193 fuel assemblies, and each 17 17 assembly contains 264 fuel rods, with three different enrichments. Both the measured data of hot zero power (HZP) and hot full power (HFP) are given in BEAVRS. The configuration of BEAVRS core is shown in Figs. 11and 12. And the parameters of core and assembly are listed in Tables 1 and 2. The detailed geometry of BEAVRS full core was built up by RMC, as shown in Figs. 13and 14. To consider the axial distributions of fission power and coolant density, the active core is divided into 10 axial segments. For thermal-hydraulics model in COBRA, only an octant of core has been considered for the 193 fuel assemblies in the full core. Therefore, the grid of 31 channels, each containing an individual fuel assembly, making up the lower triangular region of the quarter-core shown in Fig. 15. Each assembly is divided into 10 axial segments.
The HZP condition of BEARS was firstly calculated to verify the neutronics model and on-the-fly temperature treatment. The onthe-fly treatment was compared with the pre-generated cross section (XS) of exact temperature. For HZP, the coolant temperature is 566 K, the temperatures of other material (including fuel) are 293.15 K, and the boron concentration is 975 pcm. When using the pre-generated XS, XS of 293.15 K and 566 K should be prepared in advance, while only cross sections of 0 K are needed for on-thefly treatment. The results are compared in Table 3. It can be found that the Keff of on-the-fly method agrees well with the pregenerated XS. Keff of RMC also agrees well with the result of OpenMC. The relative radial pin power distributions were also compared with the results of MC21 (Kelly et al., 2013). The pin power distributions also have good agreement as shown in Fig. 16. For pre-generated XS, 200,000 neutrons per cycle are used for totally 400 cycles with 200 inactive cycles, while for on-the-fly method, 100,000 neutrons per cycle are used for totally 400 cycles with 200 inactive cycles. The ENDF/B7.0 cross-section library was used.
3.4. Results of HFP
Fig. 10. Calculations time of PWR core case.
The HFP condition of BEARS at zero burnup was calculated to verify the coupling between neutronics/thermal-hydraulics. For the coupling in HFP conditions, three feedback should be considered, including the temperature of coolant and fuel, the density of coolant and the boron concentration in coolant. For the boron concentration, there are two kinds of critical boron curves, one is in 100% core power and the other is with realistic power history, as is shown in Fig. 17. In the burnup calculations in HFP condition, the boron concentration was updated at every burnup step according to the given boron concentration of benchmark, to see whether the reactor is critical or not. However, for the HFP condition at zero burnup, there is no data of critical boron concentration, as the first burnup point in Fig. 17 for 100% power is for the fourth EFPD which is 599 pcm. Consider the effect of xenon poisoning at the beginning of burnup, the critical boron concentration must be larger than 599 pcm. Therefore, HFP condition of BEARS at zero burnup with different boron con-
440
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Fig. 12. Axial cross section BEAVRS core.
Table 1 Nominal value of thermal hydraulics condition. Item
Value
Units
Total power Initial mass flow rate Reference pressure Inlet water temperature Outlet water temperature
3411 17,083 155.17 292.78 310
MW kg/s bar °C °C
Table 2 Assembly size information. Item
Size
Number of fuel rods Number of guide tubes rods Active length (mm) Bundle pitch (mm) Fuel rod diameter (mm) Cladding inner diameter (mm) Cladding outer diameter (mm) Pin pitch (mm) Inner diameter of guide tube (mm) Outer diameter of guide tube (mm)
264 25 3657.6 215.04 7.84 8.00 9.14 12.60 11.22 12.04
Fig. 13. Radial cross section of BEAVRS core (RMC).
centrations were compared in Table 4. For HFP, 10,000 neutrons per cycle are used for totally 400 cycles with 150 inactive cycles. From Table 4, it can be found that with 975 pcm boron concentration, Keff of HZP is almost critical, while Keff of HFP is 1477.1 pcm smaller than 1. It shows the total effect of temperature feedback (including fuel, coolant and boron) is negative. On the other hand, 599 pcm is also not appropriate for HFP, as Keff is 3189.5 pcm larger. The reason is there is no xenon at zero burnup, while at the fourth EFPD the xenon poisoning reach saturation. By critical search of boron concentration, the critical boron concentration at zero burn is 850 pcm (see Fig. 18). In order to achieve the detailed power, 500,000 neutrons per cycle are used for totally 400 cycles with 150 inactive cycles. The convergence of different cycle neutrons number are compared in
Fig. 10. The ‘‘Max relative diff” (or ‘‘Max diff”) means the maximum relative difference of assembly power between two coupling iterations. It can be found that the ‘‘Max relative diff” of 500,000 neutrons is smaller. The statistic relative deviation of Monte Carlo calculations is also compared in Table 5. ‘‘Max RE” means the maximum relative deviation of assembly power in Monte Carlo calculations. It can be found that the convergence of coupling iterations depend on the statistic relative deviation of Monte Carlo calculations. When convergence is achieved, ‘‘Max diff” is within or close to three times of ‘‘Max RE”. The linear power per assembly, fuel average temperature, coolant temperature and coolant density at the middle of axial core are
441
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Fig. 17. Critical boron concentration evolution with burnup (reference). Fig. 14. Axial cross section of BEAVRS core (RMC). Table 4 Keff with different concentration.
Fig. 15. Radial channels of BEAVRS core (COBRA).
Table 3 Keff comparisons of pre-generated cross sections and on-the-fly method.
Pre-generated XS On-the-fly XS OpenMC
Keff
Diff
0.998207 (0.000103) 0.998455 (0.000178) 0.99920 (0.00004)
– 0.000248 0.000993
power
Condition
Keff
HZP + 975 ppm HFP + 975 ppm HFP + 599 ppm HFP + 850 ppm
0.998461 0.985229 1.031895 0.999847
and
boron
(0.000070) (0.000498) (0.000445) (0.000457)
shown in Figs. 19–22. It can be found that the fuel average temperature and coolant temperature have the checkerboard distributions similar to linear power distributions, proving that the results of thermal-hydraulics calculations are reasonable. The axial distributions of fuel average temperature, coolant temperature and coolant density in the hottest assembly are shown in Fig. 23. The detailed radial pin power distributions of RMC in HFP condition at zero burnup are compared with that of HZP and the results of MC21 (Daniel et al., 2014), as shown in Fig. 24. Compared with the HZP results, the power distributions of HFP are more uniform, due to the negative temperature feedback. The power distributions of RMC agree well with the results of MC21 although there are some differences. The differences are because the burnup of
RMC
MC21 Fig. 16. Relative radial pin power distributions (HZP).
442
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Another important phenomenon should be noted is that unnatural and unreasonable cross shaped power peaks appear in some hot assemblies in the middle of core. The reason is that in the thermal-hydraulics model of MC21, each assembly is divided into 4 separate channels, which gives raise to the power peaks between channels in the same assembly. 3.5. Results of burnup calculation in HFP
Fig. 18. Convergence curve of different cycle neutrons number.
Table 5 Convergence status of different cycle neutrons number. Cycle neutrons number
Max RE
Max diff
Diff/RE
500,000 10,000
0.0231 0.152
0.03139 0.26459
1.359 1.741
MC21 results is 92 calendar days (1023.3 MWD/THM) with 76% full power and 702.5 ppm boron concentration.
On the basic of neutronics/thermal-hydraulics coupling in HFP condition, the domain decomposition strategy was added to perform the full core ‘‘transport-burnup-thermal-hydraulics” of BEAVRS benchmark. Each fuel pin and burnable poison pin is treated as individual burnup regions, and each pin was divided into 10 axial segments, summing up to 640,000 burnup regions. The whole core is divided into 16 domains: 2 axial domains and 8 radial domains, as shown in Fig. 25. The large scale parallelism was performed on Tianhe2 super computer. 30 nodes with 720 processes were used, with 45 processes in each domain. For the power normalization of depletion calculation, the ratio of real power and relative power tallied by Monte Carlo code should be figured out. For domain decomposition, each domain has its total power, and communications between domains will be conducted with MPI_Allreduce to sum up the domain power to get the total tallied power. Once the total tallied power was attained, the ratio of real power and relative power can be figured out. It should be noted that in this MPI communications, only the domain power should be communicated, therefore the data size is small.
Fig. 19. Linear power per assembly (w/m).
Fig. 20. Fuel average temperature (K).
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
443
Fig. 21. Coolant temperature (K).
Fig. 22. Coolant density (g/cm3).
Fig. 23. Axial distributions.
The ‘‘transport-burnup-thermal-hydraulics” strategy was shown in Table 6. Two sequential transport/thermal-hydraulics coupling iterations were performed firstly, to attain the power distributions with thermal-hydraulics feedback considered, as shown in Fig. 26. Then the burnup calculation were performed, using the predictor-corrector depletion strategy. As the boron concentration in coolant changes in different burnup steps, the capacity of varying materials was used to change the boron concentration in each burnup step according to Fig. 17 to keep the reactor critical.
For the hot full power calculations, the specific power is 41.699 W/gHM for each burnup step. The burnup steps are: 4, 7, 5, 6, 9, 5, 16, 17, 16, 11, 14, 14 (days). . ., which is set according the reference of critical boron concentration in Fig. 16. This setting of burnup steps is convenient for the comparisons of reactivity with critical boron concentrations. To investigate the effects of thermal-hydraulics feedback in burnup calculations, another burnup calculation without thermalhydraulics feedback was performed, while changing the boron concentration according to Fig. 17. For the case without thermalhydraulics feedback, the fixed temperature cross sections were used, in which the coolant temperature is set as 600 K, and fuel temperature is 900 K. The results are compared in Fig. 27. It was found that keff of case without feedback keeps increasing when burnp increases, raising to 3% larger than critical. While keff with feedback always maintains within 0.99–1.0. It can be concluded that the critical boron concentration with thermal-hydraulics feedback by RMC agrees with the critical boron concentration given by the benchmark (Fig. 17). The difference of reactivity may be resulted from a lot of reasons. But the main reason is because of the calculations were carried on by hot full (100%) power. However, in the cycle 1 of BEAVRS, the real power is not constant but changed with burnup and the average power is about 75%. When increasing the power from 75% to 100%, keff will decrease because the temperature is higher and negative reactivity feedback. Moreover, when the power increase, the saturated xenon concentration will increase, also leading to the decrease of keff. The Relative radial pin power distributions with different burnup were calculated by RMC and compared in Fig. 28. It can be
444
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
HZP
HFP(RMC, day 0)
HFP(MC21, day 92)
Fig. 24. Relative radial pin power distributions (HFP).
Fig. 25. Diagram of domain decomposition of BEAVRS.
Table 6 Coupling strategy. Iteration
Contents
Calculation conditions
1
1st transport/1st thermalhydraulics 2nd transport /2nd thermalhydraulics Predictor depletion Corrector transport/corrector depletion
106 neutrons per cycle, 150 inactive, 100 active 106 neutrons per cycle, 150 inactive, 400 active
2 3 4
106 neutrons per cycle, 150 inactive, 400 active
found that the power distributions of core becomes uniform when the burnp increases. 4. Conclusions The‘‘transport-burnup-thermal-hydraulics” multi-physics coupling is crucial for realistic reactors simulations and benchmarks
calculations such as MIT BEAVRS benchmark. In order to perform the high fidelity‘‘transport-burnup-thermal-hydraulics” coupling, several advanced techniques were developed in RMC. The onthe-fly temperature-dependent cross section treatment was developed to consider the temperature effects on nuclear cross sections in resolved resonance region and thermal energy region. A hybrid coupling method was proposed and RMC was coupled with subchannel code COBRA to consider the distributed temperatures and coolant densities in the whole core. The hybrid coupling method can reduce the difficulty of modeling and improve the versatility of coupling by managing all the useful data by internal memory in neutronics code, while making good use of the existing thermal-hydraulics codes. The spatial domain decomposition method based on nested MPI parallelism was developed to handle the memory problem of pin-wise three-dimensional full core burnup calculations. The ‘‘transport-burnup-thermal-hydraulics” coupling system was applied to burnup simulations of BEAVRS benchmark in hot full power condition.
445
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
Fig. 27. Keff evolution with burnup.
Fig. 26. Flowchart for the transport-T/H-burnup coupling scheme.
The proposed domain decomposition method based on nested MPI parallelism was firstly verified by solving the burnup problems of a typical PWR assembly and a two-dimensional PWR core. The results of these two cases show that domain decomposition can achieve consistent kinf/keff compared with original version of RMC, while enlarging the computational burnup regions. The time ratio of with and without domain decomposition increases with burnup. The on-the-fly temperature treatment was then verified by the hot zero power condition of BEAVRS. Keff of on-the-fly method agrees well with the pre-generated cross sections of exact temperature. The relative radial pin power distributions also agree with the results of MC21.
Burnup=0
For the HFP calculations, the detailed radial pin power distributions of RMC in HFP condition at zero burnup are compared with those of HZP and the results of MC21. The critical boron concentration of HFP at zero burnup was found to be 850pcm. The convergence of coupling iterations are compared with statistic relative deviation of Monte Carlo calculations, showing that the convergence of coupling iterations depends on the statistic relative deviation of Monte Carlo calculations. The power distributions of HFP are more uniform than those of HZP due to the negative temperature feedback. Moreover, unreasonable cross shaped power peaks are found in some hot assemblies of MC21 results, which are resulted from the thermal-hydraulics model of MC21. Finally, the burnup calculations in HFP were performed. The results show that keff with feedback always maintains within 0.99–1.0 when burnup increases, while keff of case without feedback keeps increasing when burnup increases, rising to 3% larger than criticality. It can be concluded that the critical boron concentration with thermalhydraulics feedback by RMC agrees with the critical boron concentration given by the benchmark. The relative radial pin power distributions with different burnup were also compared and found that the power distributions of core become uniform when the burnup increases. This work provides the foundation for ‘‘transport-burnup-ther mal-hydraulics” coupling calculations in multi-physics coupling and lifecycle simulations of nuclear reactors using RMC code. More
Burnup=1501 MWD/THM Fig. 28. Relative radial pin power distributions with different burnup.
Burnup=3544 MWD/THM
446
S. Liu et al. / Annals of Nuclear Energy 101 (2017) 434–446
research will be done to fulfill the two cycle simulations of BEAVRS benchmark with RMC code, and more advanced thermal-hydraulic solver such as COBRA-TF would be used for coupling with pin-bypin sub-channel and parallelism in TH calculations. Acknowledgments This work is partially supported by National Science and Technology Major Project of Large Advanced PWR Nuclear Power Plant in China (2011ZX06004-024) and Project 91126017/11475098 by National Natural Science Foundation of China. References Cullen, D.E., Weisbin, C.R., 1976. Exact Doppler broadening of tabulated cross sections. Nucl. Sci. Eng. 60, 199. Daniel, K., Brian, A., Paul, R., 2014. Analysis of select BEAVRS PWR benchmark cycle 1 results using MC21 and OpenMC. In: PHYSOR2014, Kyoto, Japan, September 28–October 3, 2014. Forget, B., Xu, S., Smith, K., 2014. Direct Doppler broadening in Monte Carlo simulations using the multipole representation. Ann. Nucl. Energy 64, 78–85. Horelik, N., Herman, B., 2012. Benchmark for evaluation and validation of reactor simulations. MIT Computational Reactor Physics Group. URL:
Liang, J., Wang, K., Yu, G., et al., 2014. Research on tally data decomposition algorithms based on reactor Monte Carlo code RMC. Nucl. Power Eng. 35 (4), 142 (in Chinese). Liang, J., Wang, K., Qiu, Y., Chai, X., Qiang, S., 2016. Domain decomposition strategy for pin-wise full-core Monte Carlo depletion calculation with the reactor Monte Carlo code. Nucl. Eng. Technol. 48, 635–641. Liu, S.C., Yu, J.K., Liang, J.G., Wang, K., 2015. Study of neutronics and thermalhydraulics coupling with RMC COBRA-EN system. In: 7th International Conference on Modelling and Simulation in Nuclear Science and Engineering (7ICMSNSE), Ottawa, Ontario, Canada, October 18–21, 2015. Liu, S.C., Yuan, Y., Yu, J.K., Wang, K., 2016a. Development of on-the-fly temperaturedependent cross-sections treatment in RMC code. Ann. Nucl. Energy 94, 144– 149. Liu, S.C., Yuan, Y., Yu, J.K., Wang, K., 2016b. Reaction rate tally and depletion calculation with on-the-fly temperature treatment. Ann. Nucl. Energy 92, 277– 283. Pavlou, T.A., Ji, W., 2014. On-the-fly sampling of temperature-dependent thermal neutron scattering data for Monte Carlo simulations. Ann. Nucl. Energy 71, 411– 426. She, D. et al., 2013. Development of burnup methods and capabilities in Monte Carlo code RMC. Ann. Nucl. Energy 51, 289–294. She, D. et al., 2014. 2D full-core Monte Carlo pin-by-pin burnup calculations with the RMC code. Ann. Nucl. Energy 64, 201–205. Viitanen, T., Leppänen, J., 2012. Explicit treatment of thermal motion in continuousenergy Monte Carlo tracking routines. Nucl. Sci. Eng. 171, 165. Wang, K. et al., 2015. RMC – a Monte Carlo code for reactor core analysis. Ann. Nucl. Energy 82, 121–129. Yang, F., Liang, J.G., Yu, G.L., Wang, K., 2015. Research on on-the-fly Doppler broadening based on the stochastic algorithm. In: 2015 ANS Annual Meeting. San Antonio, TX. June 7–11. Yesilyurt, G., Martin, W.R., Brown, F., 2012. On-the-fly Doppler broadening for Monte Carlo codes. Nucl. Sci. Eng. 171, 239.