Study on the calculational framework development of the advanced numerical reactor neutronics code SHARK

The SHARK program (Simulation-based High-fidelity Advanced Reactor physics Kit) is a high-fidelity heterogeneous neutronics code for the numerical reactor system being developed at the Nuclear Power Institute of China (NPIC). The program uses a Constructive Solid Geometry (CSG) framework to model various complex geometries. To enhance the flexibility and robustness during continuous development process, SHARK program attempts to support a rich set of methods, tools and library options within a unified general framework as a “toolkit”. For the multi-core clustered supercomputer architectures that are commonly used today, the SHARK program adopts a hybrid parallel strategy of MPI and OpenMP to achieve complementary advantages between them. In addition, the framework of SHARK program is designed with a true object-oriented manner. Through reasonable abstraction, inheritance and encapsulation, the maintainability and extensibility of the code are improved, and long-term team development is facilitated. Up to now, key modules for cross-section generation, heterogeneous transport calculation and microscopic depletion have been developed under the general frameworks. The main features of SHARK’s “resonance-transport-depletion” coupling system are elaborated in this paper, and some verification and validation (V&V) results in the current phase are presented and discussed.

1 Introduction

The innovative numerical reactor concept has been proposed for several years and is continuously developed and refined. Today, numerical reactors in a general sense are often considered as a combination of key elements such as high performance computing platforms, high fidelity and high resolution modeling and simulation, multi-physics and multi-scale coupled analysis. As the core component for the digital twin of nuclear reactors, numerical reactor will play a key role in improving design, exploiting margins, exploring mechanisms, and reducing costs. It will have far-reaching influence on both existing reactors and newly developed reactors.

For numerical reactor systems with multi-physics coupling simulation, high-fidelity neutronics programs for simulating neutron/photon field distribution in the reactor are undoubtedly one of the core engines. Depending on the underlying principles, high-fidelity neutronics computational techniques can be simply divided into Monte Carlo route and deterministic theory route. In recent years, with the rising of numerical reactor, a number of excellent high-fidelity neutronics calculation programs have emerged internationally for application in this field, such as MCNP (Cardoni, 2011), MC21 (Aviles et al., 2016), RMC (Liu et al., 2018), and JMCT (Deng et al., 2018) based on Monte Carlo theory, and MPACT (Kochunas et al., 2013), NECP-X (Chen et al., 2018), DeCART (Joo et al., 2014), nTRACER (Jung, 2013), STREAM (Choi et al., 2019) based on deterministic theory. Excellent work has been done by researchers around the development and application of all these programs.

At the same time, today’s high-fidelity neutronics codes still face various challenges to some degrees. First, different types of reactors often have different fuel and structural designs, core arrangements, and control/operation modes. This requires neutronics codes to be adaptable. In terms of geometry modeling, it is worthwhile to investigate how to model complex geometries easily and intuitively, and how to partition and map different levels of meshes efficiently. Sometimes, this can also be solved partly with the help of various advanced multi-physics coupling platforms, such as SALOME, and MOOSE. In terms of fuel simulations, different designs including enrichments and dispersion manners occurs in conventional type fuel and accident resistant fuel (ATF), which can pose a variety of problems on numerical simulations. In terms of the reactor core as a whole, different spectrum characteristic, flux levels and gradients, and dynamic properties place higher demands on the cross-section library and resonance data processing, the fuel depletion chain and algorithms, and also require a more stable and powerful transport solver to support them.

Second, high-fidelity neutron transport calculation is always an extremely memory and computation intensive scientific computing task. Therefore, it is necessary to take full advantage of today’s rapidly developing high-performance computing platforms. The software parallelism strategies and their implementation techniques should be improved. On the one hand, the parallelizability of specific neutronics algorithms needs to be analyzed in depth. On the other hand, the architecture characteristics of the High Performance Computing (HPC) platform should also be fully considered.

In addition, from the point of view on software engineering, maintainability, scalability, and teamwork development are also key considerations in the software framework design and code development process for such a large-scale scientific computing software.

The SHARK program (Simulation-based High-fidelity Advanced Reactor physics Kit) is a high-fidelity heterogeneous neutronics code for numerical reactor systems being developed at the Nuclear Power Institute of China (NPIC). Like other similar programs, the team is trying to achieve a step-by-step exploration and solution of the aforementioned problems during the development. It uses the popular Constructive Solid Geometry (CSG) strategy for geometric modeling, and an ordered hierarchical structure for grid management and mapping between modules of neutronics calculations. Thus, a wide range of complex problem geometries can be described through text cards or visualization tools under development. The SHARK program attempts to enhance the adaptability and robustness as a “toolkit” by supporting a rich and flexible set of methods, tools and library options within a unified overall framework. For the multi-core clustered supercomputer architectures that are commonly used today, it adopts a hybrid parallelism strategy of MPI and OpenMP to achieve complementary advantages. In addition, the framework of SHARK program is designed with a true object-oriented philosophy. Through reasonable abstraction, inheritance and encapsulation, it improves the maintainability and extensibility of the code and facilitates long-term collaborative team development.

SHARK program is currently being developed continuously according to the above concepts. Up to now, the parts that have been developed or integrated include the main framework, the cross-section library and resonance module, the heterogeneous transport module, the high-fidelity depletion module and relevant utilities. The numerical simulation capability of steady-state and dynamic depletion for heterogeneous neutron transport problems at the microscopic cross-section level has been basically established. For the treatment of resonance self-shielding effects, the SHARK framework currently offers the option of improved subgroup methods (Zhang et al., 2022a) with corresponding cross-section library, and other options like global-local method (He et al., 2018) are also being integrated. For heterogeneous transport calculations, the SHARK program uses the method of characteristics (MOC) and prepares two-dimensional/one-dimensional (2D/1D), quasi-three-dimensional (Quasi-3D) (Zhao et al., 2021) and three-dimensional (3D) options. For the high-fidelity simulation of burnup, the SHARK program supports both power and pure decay modes and considers the burnup characteristics of poisons with rapid reaction rate changes. While using a refined burnup chain, the SHARK program is also designed to support easy switching of the depletion chain and library in order to meet the potential needs of different application scenarios. In addition, the development of functions such as multi-physics coupling, transient calculation, engineering calculation sequences, and output visualization are also underway or forthcoming.

This paper presents the establishment of the “resonance-transport-depletion” coupling calculation system of SHARK, and gives some V&V results of the current phase. The remaining of this paper can be divided into three parts. In Section 2, we present the general framework of the SHARK program, the methodology and technical solutions of the main computational modules. In Section 3, we give some staging results and data analysis of the program. In the last section, we conclude the paper and give some discussions.

2 Methodologies

2.1 General framework

CSG strategies are well known in computer industrial design. In the field of neutronics computing, CSG was first applied in various Monte Carlo software and has recently been applied in deterministic neutron transport software (Boyd et al., 2014). The SHARK program uses the CSG strategy to describe complex geometric models including two-dimensional and three-dimensional problems. In the CSG model, various surfaces are the most basic constructive elements and they are represented by equations. The geometric modeling module currently used by SHARK supports a variety of planes and quadratic surfaces. By defining surfaces, a number of halfspaces can be partitioned. By performing Boolean operations such as intersection, union, and complement to the halfspaces, the basic unit cell of the CSG geometry is obtained. The entity cells can be filled with material or geometric units. A large number of solid cells form universes by Boolean operations. In addition, repeated structures in the reactor can be easily and quickly defined by means of lattices. In conclusion, hierarchical geometry framework of CSG strategy is highly universal and extensible, which can effectively enhance the front-end geometry modeling capability of the program. At this stage, the SHARK program has successfully tried to model reactor geometry in the form of square and hexagonal lattices with fuel rods or plates (Figure 1).

FIGURE 1

FIGURE 1. CSG modeling function in SHARK program.

The logic flow of the high-fidelity neutronics calculation program is relatively complex and has a large number of program modules and coupling interfaces. Therefore, in order to facilitate the integration of different models within the SHARK program framework, we strictly follow the object-oriented programming ideas. In the process of building the main calculational framework, the basic principle of “high cohesion - low coupling” is used to cut the key modules. At the same time, a large number of utility modules are provided to help program developers reuse them. For the different method implementations in the same computational stage, the concept of inheritance between parent and child classes is implemented. For example, “Resonance” class is the parent class of resonance self-shielding computation, which defines the abstraction process, basic data and methods in resonance calculations. “SubGroup” class, “Global-Local” class and other specific resonance methods are the child classes of “Resonance”. They inherit the basic features of the “Resonance” class, but also define different implementation details in an individualized way by static or dynamic polymorphism. The different child classes can be developed, extended and maintained independently under a unified framework. Figure 2 shows a simple illustration of this design philosophy, where a mix of C++ and Python programming is used in the development of SHARK program. In particular, C++ is used to develop the calculational kernel, which facilitates stable and efficient numerical computation. The application interface, on the other hand, is written in Python, which is more user-friendly and extensible.

FIGURE 2

FIGURE 2. Design of the SHARK program framework for depletion.

MPI and OpenMP hybrid parallel strategies are known to perform well in today’s commonly adopted distributed-shared memory supercomputer platforms. MPI technology is used to handle process-level parallelism for distributed memory with good scalability. OpenMP technology is used to handle thread-level parallelism for shared memory with simple implementation and almost without communication latency and load balancing issues. On high performance computing, the SHARK program provides this hybrid parallel architecture. Specifically, the MPI technique is used to handle spatial domain decompositions, with each domain independently performing ray tracing and communicating with neighboring domains on the inner boundary (Figure 3); the OpenMP technique is used to parallelize the sweep of rays, meshes and energy groups. This hybrid parallel strategy can effectively distribute the computation and storage load with good parallel efficiency.

FIGURE 3

FIGURE 3. Domain decomposition parallel strategy in SHARK program.

In SHARK program, data is carefully classified and encapsulated. Data inside a class is only available for manipulation by trusted classes. Various pointers are also used in the program implementation to improve data manipulation efficiency and save memory cost. In addition, for some special data, linked lists, binary trees, hash tables and the corresponding sorting and searching algorithms are used effectively.

To ensure compatibility, the programming language standard is used as the basic programming principles. This minimizes the dependence on specific compilers and compilation environments. Some settings that may be related to the system environment are handled using conditional compilation. In addition, common tools for version control and automatic testing are also essential in the code development.

2.2 Cross-Section generation

Cross-section generation plays an important role in neutronics calculations. Especially for a program like SHARK, which is based on the one-step heterogeneous transport theory, the largest source of error and uncertainty in the system comes from the multi-group cross-section generation. It directly affects the accuracy and applicability of the program. The cross-section generation module mainly consists of two parts: resonance calculation and the corresponding nuclear data library. In addition to accuracy and efficiency, the SHARK program imposes other requirements on the resonance effect treatments. On the one hand, it has to be able to provide spatially dependent high-resolution results for multi-physics coupling within the fuel element. On the other hand, it has to be highly geometrically adaptable to meet the application requirements of different reactor types. The subgroup resonance option inside the SHARK program is illustrated (Zhang et al., 2022a).

In the subgroup approach, the effective resonance cross section ${\sigma }_x$ is expressed in the form of

${\sigma }_x=\frac{{\displaystyle \sum }_n{w}_{x,n}{\sigma }_{x,n}{\varphi }_n}{{\displaystyle \sum }_n{w}_{a,n}{\varphi }_n}(1)$

where $\left(w_{x,n},{\sigma }_{x,n}\right)$ is the x-type subgroup parameter; ${\varphi }_n$ is the subgroup flux.

In Eq. 1, if it is required to take ${\sigma }_x$, the subgroup equations with the form of transport equations must be solved to obtain the subgroup fluxes. To improve the efficiency of the core-level resonance calculation, the approximate solution (Stamm’ler, 2001) of the slowing down equation is introduced to rewrite Eq. 1 as

${\sigma }_x\approx \frac{{\displaystyle \sum }_n{w}_{x,n}{\sigma }_{x,n}\frac{{\sigma }_{b,n}}{{\sigma }_{a,n}+{\sigma }_{b,n}}}{{\displaystyle \sum }_n{w}_{a,n}\frac{{\sigma }_{b,n}}{{\sigma }_{a,n}+{\sigma }_{b,n}}}(2)$

${\sigma }_{b,n}=\frac{\lambda {\Sigma }_p+{\Sigma }_{e,n}}{N}(3)$

${\Sigma }_{e,n}={\Sigma }_{b,n}-\lambda {\Sigma }_p=N\left(\frac{{\varphi }_n}{1-{\varphi }_n}\right){\sigma }_{a,n}-\lambda {\Sigma }_p(4)$

where ${\sigma }_{b,n}$ is the subgroup background cross section; $\lambda {\Sigma }_p$ is the product of the intermediate resonance factor and the macroscopic potential scattering cross section; ${\Sigma }_e$ is the macroscopic equivalent cross section; and $N$ is the resonant nuclide number density.

The key to calculate the effective resonance cross section using Eqs. 2–4 is to obtain ${\Sigma }_{e,n}$, which is also uneconomical to solve directly. Since ${\Sigma }_e$ and ${\sigma }_a$ present a relatively smooth monotonic relationship, a linear logarithmic interpolation table of ${\Sigma }_e\left({\sigma }_a\right)$ can be constructed to reduce the frequency of solving the subgroup equations.

The interpolation table of subgroup equivalent cross sections ${\Sigma }_e\left({\sigma }_a\right)$ is related to resonance energy groups and resonant nuclides. According to the distribution of the resonance peaks, two strategies, " resonant nuclide classification” and “resonance energy segmentation”, can be used to reduce the computational cost.

Based on the above simplification, the subgroup equations as Eq. 5 are solved by MOC to obtain the subgroup fluxes and complete the tabulation.

$\begin{array}{cc}\mathbf{\Omega }\cdot \nabla {\psi }_m\left(\mathbf{r},\mathbf{\Omega }\right)+\left({\Sigma }_{a,m}+{\Sigma }_p\right){\psi }_m\left(\mathbf{r},\mathbf{\Omega }\right)=\lambda {\Sigma }_p+\left(1-\lambda \right){\Sigma }_p{\varphi }_m\left(\mathbf{r}\right),& m=1,...,M\end{array}(5)$

where, $\mathbf{\Omega }$ is the direction vector; $\mathbf{r}$ is the position vector; ${\psi }_m$ is the subgroup m angular flux.

When multiple resonant nuclides are contained in the same material region, the iterative idea is used to deal with the multi-nuclide resonance interference effect.

The effective resonance absorption integral after convergence of the Bondarenko iteration can be used to obtain the background cross section by the binary search method, and the effective resonance scattering cross section is calculated accordingly.

For non-resonant nuclides, the slowing down effect of light nucleus is relatively accurately treated by correction factor to take into account resonance spectrum effects on the scattering cross section of claddings and moderators.

The subgroup parameters and multigroup cross sections corresponding to the subgroup module are provided by the accompanying cross section library (Chai et al., 2017). At present, the library contains more than 280 nuclides in total, considering the main reaction types such as fission, capture, (n, 2n)/(n, 3n), and decay. The number of its energy groups will also cover a varying range to meet the application needs of different scenarios.

2.3 Transport calculation

Using the generated multigroup cross sections as input, the transport solver completes the calculation and provides the flux distributions for a specific condition. The SHARK program uses MOC with powerful heterogeneous neutron transport capability as the solver. Besides direct 3D MOC, 2D/1D MOC and Quasi-3D MOC (Zhao et al., 2021) have been developed. In general, both of them transform the 3D transport problem into a 2D MOC problem with several layers coupled to each other in the axial direction.

In 2D/1D MOC, the 3D neutron transport equations are integrated axially and radially to obtain the basic coupling equation as

${\xi }_m\frac{\partial {\psi }_{g,m}^L\left(x,y\right)}{\partial x}+{\eta }_m\frac{\partial {\psi }_{g,m}^L\left(x,y\right)}{\partial y}+{\Sigma }_{t,g}^L\left(x,y\right){\psi }_{g,m}^L\left(x,y\right)=Q_{g,m}^L\left(x,y\right)-T{L}_{g,m,L}^{Axial}\left(x,y\right)(6)$

${\mu }_m\frac{d{\psi }_{g,m}^p\left(z\right)}{dz}+{\Sigma }_{t,g}^p\left(z\right){\psi }_{g,m}^p\left(z\right)=Q_{g,m}^p\left(z\right)-T{L}_{g,m,p}^{Radial}\left(z\right)(7)$

where, ${\psi }_{g,m}^L\left(x,y\right)$ is the angular flux of the direction m energy group g at the position of $\left(x,y\right)$ on the axial Lth layer; ${\psi }_{g,m}^p\left(z\right)$ is the angular flux of the direction m energy group g in the radial mesh p at the axial position of z; $T{L}_{g,m,L}^{Axial}\left(x,y\right)$ and $T{L}_{g,m,p}^{Radial}\left(z\right)$ is the axial leakage term and radial leakage term, respectively; the rest of the symbols are in common notation.

The SHARK program uses a discrete-ordinates method (S_n) in finite difference format to solve the 1D Eq. 7 in order to fully consider the heterogeneous transport effects in the axial direction.

For the potential numerical instability problem in 2D/1D MOC calculation, the isotropic leakage approximation and leakage splitting methods (Zhao et al., 2018) are used.

And in the Quasi-3D MOC, the 3D neutron transport equation is transformed into the 2D coupling equation by introducing the axial differential format approximation as follows.

$\begin{array}{c}\epsilon \frac{\partial {\psi }_{g,m}^L\left(x,y\right)}{\partial x}+\eta \frac{\partial {\psi }_{g,m}^L\left(x,y\right)}{\partial y}+\left({\Sigma }_{t,g}^L\left(x,y\right)+\frac{2\mu }{\Delta z_L}\right){\psi }_{g,m}^L\left(x,y\right)=Q_{g,m}^L\left(x,y\right)+\frac{2\mu }{\Delta z_L}{\psi }_{g,m}^{L-\mathit{1}/\mathit{2}}\left(x,y\right)\mu >0\\ \epsilon \frac{\partial {\psi }_{g,m}^L\left(x,y\right)}{\partial x}+\eta \frac{\partial {\psi }_{g,m}^L\left(x,y\right)}{\partial y}+\left({\Sigma }_{t,g}^L\left(x,y\right)-\frac{2\mu }{\Delta z_L}\right){\psi }_{g,m}^L\left(x,y\right)=Q_{g,m}^L\left(x,y\right)-\frac{2\mu }{\Delta z_L}{\psi }_{g,m}^{L+\mathit{1}/\mathit{2}}\left(x,y\right)\mu <0\end{array}(8)$

${\psi }_{g,m}^L\left(x,y\right)=\frac{{\psi }_{g,m}^{L+\mathit{1}/\mathit{2}}\left(x,y\right)+{\psi }_{g,m}^{L-\mathit{1}/\mathit{2}}\left(x,y\right)}{2}(9)$

where ${\psi }_{g,m}^{L\pm \mathit{1}/\mathit{2}}\left(x,y\right)$ is the angular flux of the direction m energy group g at the positions $\left(x,y\right)$ on the top and bottom surfaces of the Lth layer; the rest of the symbols are in common notation.

As seen in Eqs. 8, 9, the Quasi-3D MOC does not contain a transverse leakage term and the right term of the equation is always positive, thus providing better numerical stability than the 2D/1D MOC. All these transport solvers are accelerated by the Coarse Mesh Finite Difference (CMFD) method.

2.4 Depletions

The depletion calculation is used to track the nuclides over time and predict the dynamic behavior of the core during burnup. It is important for fuel management, breeding and transmutation, spent fuel reprocessing, and multi-physics coupling analysis of the core. In the framework design of the high-fidelity transport software, the depletion function is not only responsible for completing the burnup calculation, but also driving the solvers of resonance, transport, and multi-physics to complete the coupling analysis of each burnup point (Figure 2).

The SHARK program adopts an internally coupled framework for depletion calculations (Zhang et al., 2022b), with the following features reserved.

1) Support for both neutron reaction depletions with multiple burnup steps and pure decay calculations during the cooling down stage on the reactor shutdown.

2) Support replacement for different depletion libraries and chains. As mentioned earlier, different application scenarios in various reactor types may have different requirements for depletion accuracy and efficiency. The implementation of many codes is usually tied to a specific chain, which makes it difficult to achieve easy replacement of depletion data. The SHARK program is designed to load the depletion chain and depletion data as external library files for “plug-and-play” scalability.

3) Several “transport-depletion” coupling schemes and point depletion solver options are reserved. For common fuels, a conventional predictor-corrector (PC) with sub-steps strategy is used. For fuels containing special poisons such as Gd, a combination of PC and log-linear reaction rate (LLR) (Carpenter and WOLF, 2010) is used. Also, options such as Semi-PC plus LLR are considered for future multi-physics coupling. Solvers for point depletion equations with good performance, such as CRAM (Pusa and Leppanen, 2009) and Krylov subspace methods (Yamamoto et al., 2007), are also available to program developers and users in the form of different child classes in Bateman solver parent class.

4) With the ability of mapping transport meshes to depletion meshes of different sizes, the parallel capability based on domain decomposition to improve the computational efficiency of large-scale high-fidelity depletion is adopted.

5) Provide engineering analysis functions related to depletions, including settings for power, control rods, core inlet coolant temperatures. Functions of xenon mode, soluble boron depletion, and restart calculation are also presented.

6) The depletion matrix is processed in a proper manner to avoid redundant matrix generation and excessive storage, improving the data access efficiency. In the matrix utilities, the sparse matrix storage format and operation methods suitable for different point depletion solvers are provided.

3 Preliminary V&V results

At present, the main framework of SHARK program has been established, and the first batch of key calculation modules have been developed and embedded into it. The program has the initial capability to perform the “resonance-transport-depletion” calculation, and preliminary verifications have been carried out. Due to space limitation, some of the results are shown in this paper. More details can be found in published or forthcoming articles. All the numerical results are calculated with the option of “subgroup-2D/1D”. Accordingly, a 45-group cross-section library and out-flow transport correction (Stamm’ler and Abbate, 1983) are used.

3.1 VERA depletion benchmark problems

The series of VERA progression benchmarks is one of the most authoritative benchmark problems for validating heterogeneous transport calculation codes (Godfrey, 2014; Kang, 2015). Among them, the VERA1/2 series are for commercial Pressurized Water Reactor (PWR) pins and lattices. They consider heterogeneous effects from fuel enrichment, temperature, burnable poisons (IFBA/WABA/gadolinium rods), control rods (AIC/B₄C), and spacer grids. We have already completed the steady-state verifications in our previous study and obtained good results (Zhang et al., 2022c). In this paper, we focus on verifying the newly developed depletion functions, where the reference solution (Wen, 2022) is obtained from the advanced high-fidelity transport calculation program MPACT. For the single pin problems, the calculation scheme of 10-6-0.01 (10 azimuthal and 6 polar angles in [0,π] with 0.01 cm ray width) is used. The scheme of 16-6-0.01 is used in the single lattice problems.

The deviations of the infinite multiplication factor (k_inf) with depletion for the single pin problems are shown in Figure 4. As can be seen, all the results agree well with the reference solution, except for the 1E calculation. For different enrichments and temperatures, the overall deviation falls within the range of -300 pcm to +100 pcm and varies smoothly with the burnups. For the 1E problem (IFBA pin with 3.1w/o enrichment), the SHARK program shows a large abnormal deviation from the reference solution; however, after comparing with the Monte Carlo program Serpent (Leppänen, 2013), it is found that SHARK and Serpent agree well and the trend of deviations with burnups is normal.

FIGURE 4

FIGURE 4. Deviation of infinite multiplication factor for VERA1 depletion problem.

For the single lattice problems, the k_inf and its deviation with burnups for 2C (3.1 w/o enrichment), 2G (24 AIC control rod insertion), and 2M cases (128 IFBA) are given in Figure 5. As can be seen from the figure, for the lattices with and without absorbers, the calculation results of the SHARK program agree well with the reference solutions.

FIGURE 5

FIGURE 5. Infinite multiplication factors and deviations for the VERA2 problems.

Figure 6 shows the deviation statistics of the pin power distributions during depletion. The deviation of pin powers does not exceed 0.4% in the case without control rod, and does not exceed 0.8% in the case with inserted control rods. Moreover, the deviations are stable even in deep burnups.

FIGURE 6

FIGURE 6. Deviations of pin power distributions in the VERA2 problems.

3.2 VERA4 3-D colorset benchmark problem

The problem is a 3-D multi-assembly one for PWR in a 3 × 3 arrangement (Godfrey, 2014). It takes into account different enrichments, burnable poisons (Pyrex), different AIC/B₄C control rod insertions, and structural details such as end plugs, springs, nozzles, and spacer grids, which can truly reflect the software’s ability to simulate radial and axial heterogeneous effects (Figure 7).

FIGURE 7

FIGURE 7. Configurations of the VERA4 3D problem.

The results of SHARK are given in Table 1 with the reference from MPACT, which are obtained by the scheme of 16-6-0.03 with about 50 planes in the active zone. Both the Differential Rod Worth (DRW) and Integrated Rod Worth (IRW) are very close to the reference solutions. In general, the maximum deviation is +28pcm for DRW and -31pcm for IRW.

TABLE 1

TABLE 1. Calculation results for the VERA #4 problem.

The AO values also match well in Table 1. The axial power distributions at different insertions are given in Figure 8. It can be seen that the spacer grid effect and the axial power shift due to the insertion of control rods are well represented.

FIGURE 8

FIGURE 8. Axial power distributions at different control rod positions.

3.3 Macroscopic BEAVRS core problem

The problem is adapted from a large commercial PWR core loaded with 193 UO₂ assemblies (Zhang et al., 2022c). The core is arranged with fuel assemblies of 1.6 w/o, 2.4 w/o, and 3.1 w/o enrichments. The core has simplified geometry with an active zone of 380 cm high and 40 cm reflectors at both axial end (Figure 9).

FIGURE 9

FIGURE 9. Configurations of the BEAVRS core.

The reference solution of the problem is obtained from the famous Monte Carlo program RMC (Wang et al., 2013). In the calculations, the scheme of 16-6-0.02 with 26 planes is used. The convergence criteria of eigenvalue and source term are 1 pcm and 10-4 respectively. A summary of the results is shown in Table 2, and the pin power distribution deviations are shown in Figure 10. The simulation takes about 0.88 h with 578 CPU cores, in which the time partitions of track generation, MOC solution and CMFD solution are about 10%: 80%: 10%.

TABLE 2

TABLE 2. Calculation results of macroscopic BEAVRS core problem.

FIGURE 10

FIGURE 10. Radial pin power deviations in Macroscopic BEAVRS core problem.

3.4 Advanced PWR whole-core problem

This whole-core problem is from the third-generation million-kilowatt-level commercial PWR HuaLong (Figure 11). The active zone of this core consists of 177 AFA3G fuel assemblies with three different enrichments, and Pyrex is used as burnable absorbers. The SHARK program is used to simulate the three-dimensional core under the Hot Full Power (HFP), ARO condition. The heterogeneous pins with different compositions are modelled; the radial and axial reflectors outside the active zone are considered explicitly in the simulation. The scheme of 16-6-0.03 with 26 planes is used. There are tens of millions of meshes in total and 578 CPU cores are used. Using the Monte Carlo program RMC as the reference, the statistical results are given in Table 3. Assembly power distribution and their deviations are shown in Figure 12. As can be seen, SHARK has good calculation results for the initial state core of the real nuclear power plant. Only 88 pcm deviation exists for k_eff. The power distribution maintains symmetry and has a Max/Min deviation no more than ±2.6%. The improved subgroup resonance computation takes less than 20% of the total resonance-transport solution time.

FIGURE 11

FIGURE 11. Configurations of advanced PWR whole-core problem.

TABLE 3

TABLE 3. Calculation results of HuaLong whole-core problem.

FIGURE 12

FIGURE 12. Assembly power distribution and deviations in advanced PWR whole-core problem.

4 Conclusion

This paper presents the design and development of the high-fidelity neutronics calculation program SHARK with main technical features, and gives a series of stage verification results from pin to whole-core level. The results show that the SHARK program has the initial computational capability to obtain accurate and reasonable results for microscopic problems in the commercial PWRs. In the future work, the improvement and expansion of functional modules will be carried out continuously for better adaptability to different application targets and scenarios under the unified frameworks. In addition, further V&V based on measurements of nuclear power plants will move on.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

ZH: Framework design; resonance module, depletion module, utilities development; VERA verifications. ZW: Team leader, technical advisor. ZC: transport module development; BEAVRS verifications. WB: PWR Core verifications, utilities development. CZ: Project management, technical advisor. PX: Project management, technical advisor. LQ: Project management, technical advisor. YY: Project management, technical advisor. GZ: Project management, technical advisor. ZW: Project management, technical advisor.

Funding

This work is supported by National Natural Science Foundation of China (Grant No. 11905214, 12005214) and China Association for Science and Technology (Young Elite Scientists Sponsorship Program 2019QNRC001).

Acknowledgments

The first author would like to thank Wen Xingjian for sharing reference solutions in VERA depletion verifications.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Aviles, B., Kelly, D., and Aumiller, D. (2016). “Coupled MC21 and COBRA-IE solution to VERA core physics benchmark problem# 6[C],” in PHYSOR 2016, Sun Valley, Idaho, USA, May 1, 2016 - May 5, 2016.

Google Scholar

Boyd, W., Shaner, S., Li, L., Forget, B., and Smith, K. (2014). The OpenMOC The OpenMOC method of characteristics neutral particle transport codeethod of characteristics neutral particle transport code. AnnAnn. Nucl. Energy Nucl. Energy 68, 43-52. doi:10.1016/j.anucene.2013.12.012

CrossRef Full Text | Google Scholar

Cardoni, J. N. (2011). “Nuclear reactor multi-physics simulations with coupled MCNP5 and STAR-CCM+[D],” in M&C 2011: International conference on mathematics and computational methods applied to nuclear science and engineering, Rio de Janeiro, Brazil, 8-12 May 2011 (Illinois, USA: University of Illinois at Urbana-Champaign).

Google Scholar

Carpenter, D. C., and Wolf, J. H. (2010). “The log linear rate constant power depletion method [C],” in PHYSOR 2010, Pittsburgh, PA, 9-14 May 2010.

Google Scholar

Chai, X., Tu, X., Guo, F., Yin, Q., Shi, H., Wei, L., et al. (2017). Development and preliminary V&V for advanced neutron transport lattice code KYLIN-2 [J]. High Power Laser Part. Beams (in Chinese) 29 (01), 016016. doi:10.11884/HBLPB201729.160306

CrossRef Full Text | Google Scholar

Chen, J., Liu, Z. Y., Zhao, C., He, Q., Zu, T., Cao, L., et al. (2018). A new high-fidelity neutronics code NECP-X. Annals of Nuclear Energy 116, 417-428. doi:10.1016/j.anucene.2018.02.049

CrossRef Full Text | Google Scholar

Choi, S., Choe, J., Nguyen, K., Lee, W., Kim, W., Kim, K., et al. (2019). “Recent development status of neutron transport code STREAM,” in Transactions of the Korean Nuclear Society Spring Meeting, Jeju, Korea, May 23-24, 2019.

Google Scholar

Deng, L., Li, G., Zhang, B. Y., Shangguan, D., Hu, Z., Li, R., et al. (2018). “JMCT v2.0 Monte Carlo code with integrted nuclear system feedback for simulation of BEAVRS model[C],” in PHYSOR 2018, Cancun, Mexico, April 22-26, 2018.

Google Scholar

Godfrey, A. T. (2014). VERA core physics benchmark progression problem specifications, revision 4[R]. CASL, CASL-U-2012-0131-004.

Google Scholar

He, Q., Cao, L., and Liu, Z. (2018). Research of global-local resonance self-shielding calculation method based on NECP-X[J]. Nuclear Power Engineering 39 (2), 124–128. doi:10.13832/j.jnpe.2018.02.0124

CrossRef Full Text | Google Scholar

Joo, H. G., Cho, J., Kim, K., Lee, C., and Zee, S. (2014). “Methods and performance of a three-dimensional whole-core transport code DeCART,” in PHYSOR 2004, Chicago, Illinois, April 25-29, 2004.

Google Scholar

Jung, Y. S. (2013). “Development of practical numerical nuclear reactor for high fidelity core analysis[D],” (Korea: Seoul National University). Thesis(Dr. Eng.).

Google Scholar

Kang, S. (2015). Specification for the VERA depletion benchmark suite, revision 0[R]. CASL, CASL-X-2015-1014-000.

Google Scholar

Kochunas, B., Collins, B., Jabaay, D., Downar, T., and Martin, W. (2013). “Overview of development and design of MPACT: Michigan parallel characteristics transport code,” in M and C 2013: 2013 International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering, Sun Valley, Idaho, May 5-9,2013.

Google Scholar

Leppänen, J. (2013). Serpent – a continuous-energy Monte Carlo reactor physics burnup calculation code[R]. Espoo, Finland: VTT Technical Research Centre of Finland.

Google Scholar

Liu, S. C., Luo, Z., Guo, X. Y., Yu, G., Li, Z., and Kan, W. (2018). “RMC/CTF multiphysics solutions to VERA benchmark problem 6 and 7[C],” in PHYSOR 2018, Cancun, Mexico, April 22-26, 2018.

Google Scholar

Pusa, M., and Leppanen, J. (2009). Computing the Computing the Matrix Exponential in Burnup Calculationsatrix exponential in burn-up calculations[R]. Nucl. Sci. Eng. 164, 140. 150. doi:10.13182/NSE09-14

CrossRef Full Text | Google Scholar

Stamm'ler, R. (2001). HELIOS methods[R]. Nyköping, Sweden: Studsvik Scandpower.

Google Scholar

Stamm’ler, R. J., and Abbate, M. J. (1983). Methods of steady-state reactor physics in nuclear design. London: Academic Press.

Google Scholar

Wang, K., Li, Z., and Liang, J. (2013). Rmc – a Monte Carlo code for reactor core analysis[J]. La Grange Park (United States): American Nuclear Society - ANS.

Google Scholar

Wen, X. (2022). Numerical reactor high-fidelity cycle depletion and source term calculation method development and application[D]. Xi'An, Shaanxi, China: Xi’an Jiaotong University. In Chinese.

Google Scholar

Yamamoto, A., Tatsumi, M., and Sugimura, N. (2007). Numerical Numerical Solution of Stiff Burnup Equation with Short Half Lived Nuclides by the Krylov Subspace Methodolution of stiff bumup equations with short Half lived nuclides by the Krylov subspace method[J]. J. Nucl. Sci. Technol. 44 (2), 147. 154. doi:10.1080/18811248.2007.9711268

CrossRef Full Text | Google Scholar

Zhang, H. B., Peng, X. J., and Zhao, C. (2022). Subgroup method for the high fidelity neutronics code SHARK and preliminary benchmarking[J]. Nuclear Engineering and Design. under review.

Google Scholar

Zhang, H. B., Zhao, C., Peng, X. J., and Zhao, W. B. (2022c). Numerical reactor high-fidelity neutronics program SHARK development[J]. Atomic Energy Science and Technology 56 (2), 2022. In Chinese. doi:10.7538/yzk.2021.youxian.0902

CrossRef Full Text | Google Scholar

Zhang, H. B., Zhao, W. B., and Zhao, C. (2022). “Technical features of the depletion framework in the digital reactor neutronics code SHARK [C],” in Pacific Basin Nuclear Conference, Chengdu, China, Oct 31-Nov 4, 2022. under review.

Google Scholar

Zhao, C., Liu, Z., Liang, L., Chen, J., Cao, L., and Wu, H. (2018). Improved leakage splitting method for the 2D/1D transport calculation. Progress in Nuclear Energy 105 (5), 202-210. doi:10.1016/j.pnucene.2018.01.007

CrossRef Full Text | Google Scholar

Zhao, C., Peng, X. J., Zhang, H. B., Zhao, W., Li, Q., and Chen, Z. (2021). Analysis and comparison of the 2D/1D and quasi-3D methods with the direct transport code SHARK. Nuclear Engineering and Technology 54, 19. 29. doi:10.1016/j.net.2021.07.038

CrossRef Full Text | Google Scholar