Abstract
ARES is a multidimensional parallel discrete ordinates particle transport code with arbitrary order anisotropic scattering. It can be applied to a wide variety of radiation shielding calculations and reactor physics analysis. ARES uses state-of-the-art solution methods to obtain accurate solutions to the linear Boltzmann transport equation. A multigroup discretization is applied in energy. The code allows multiple spatial discretization schemes and solution methodologies. ARES currently provides diamond difference with or without linear-zero flux fixup, theta weighted, directional theta weighted, exponential directional weighted, and linear discontinuous finite element spatial differencing schemes. Discrete ordinates differencing in angle and spherical harmonics expansion of the scattering source are adopted. First collision source method is used to eliminate or mitigate the ray effects. Traditional source iteration and Krylov iterative method preconditioned with diffusion synthetic acceleration are applied to solve the linear system of equations. ARES uses the Koch-Baker-Alcouffe parallel sweep algorithm to obtain high parallel efficiency. Verification and validation for the ARES transport code system have been done by lots of benchmarks. In this paper, ARES solutions to the HBR-2 benchmark and C5G7 benchmarks are in excellent agreement with published results. Numerical results are presented which demonstrate the accuracy and efficiency of these methods.
1. Introduction
Particle transport problems arise in many different areas of engineering physics. There are two main types of simulation approaches in particle transport modeling: stochastic (Monte Carlo) and deterministic [1]. Deterministic radiation transport has gained popularity in recent years as a consequence of continuous advancements in computer technology and numerical algorithm.
ARES [2] is a multidimensional parallel discrete ordinates particle transport code that uses state-of-the-art solution methods to obtain accurate solutions to the Boltzmann transport equation. The ARES transport code system consists of seven main modules: DONTRAN1D, DONTRAN2D, DONTRAN3D, RAY2D, RAY3D, ARES_PRE, and ARES_POST. DONTRAN1D, DONTRAN2D, and DONTRAN3D are the series of DONTRAN that solve the one-dimensional, two-dimensional, and three-dimensional transport problems, respectively. RAY adopts first collision source method for ray effects mitigation. To provide the corresponding data for the transport calculation, we developed the preprocessing module ARES_PRE, which can dispose the geometry and material information of the calculated model and deal with the quadrature sets and cross-sectional message. In addition, it can also calculate the reactor core fixed source and provide the interface for the Monte Carlo and discrete ordinates coupled code. ARES_POST module presented the function to handle the calculated flux, including calculating the dose equivalent rate, DPA, and fast neutron fluence using the neutron spectrum adjustment method. The ARES can be applied to reactor physics, reactor pressure vessel fluence calculations, radiation shielding and protection, spent fuel storage cask design and analysis, fusion neutronics analysis, and criticality safety calculations.
Particle transport is an extremely challenging computational problem since the governing equation is six-dimensional with a high degree of coupling between these variables. A multigroup discretization is used in energy. Discrete ordinates differencing in angle and spherical harmonics expansion of the scattering source are adopted. The angular variable is usually discretized by replacing angular integrals with quadrature sums. The code allows multiple spatial discretization schemes and solution methodologies. A variety of spatial differencing scheme options are available, including diamond difference with or without linear-zero flux fixup, theta weighted, directional theta weighted, exponential directional weighted, and linear discontinuous finite element scheme. This discretization produces a large, sparse, linear system of equations in a seven-dimensional phase space.
The most general solution technique is source iteration. However, for optically thick problems dominated by scattering, the source iteration method converges very slowly. ARES uses the Koch-Baker-Alcouffe parallel sweep algorithm to obtain high parallel efficiency.
Since discrete ordinates method was put forward by Carlson, many computer codes based on have been developed during the past half century. In the 1960s, Los Alamos National Laboratory developed the first modern code DTF-IV [3], which considers only one-dimensional transport and employs diamond differencing with fixup. Following DTF-IV, TWOTRAN and THREETRAN [4], which can deal with 2D geometry and 3D Cartesian geometry, were developed. As a successor to DTF-IV, ONETRAN included 1D linear discontinuous spatial differencing in the 1970s. After 1980, codes with DSA method emerged, such as ONEDANT, TWODANT, and THREEDANT [5]. Corresponding to Los Alamos National Laboratory, Oak Ridge National Laboratory developed ANISN, DORT, and TORT [6]. Since the 1990s, the rapid developments in computer technology have made it available to implement massive parallel. PARTISN [7], a parallel 3D rectangular mesh neutral particle transport code, was developed to solve first-order transport equation. ATTILA [8], based on unstructured tetrahedral mesh, uses linear discontinuous finite element approximation in the spatial variables. In the 2000s, Oak Ridge National Laboratory developed Denovo [9] to replace TORT. Denovo, which has multiple spatial differencing schemes, uses KBA parallel sweep algorithm and DSA preconditioned Krylov solver. Recently, THOR [10] implemented Arbitrarily High Order Transport method of the Characteristic type.
This paper is organized as follows. In Section 2, the spatial, angular, time, and energy discretizations used in ARES are described and summarized. The parallel algorithms and solution techniques will be discussed in Section 3. The first collision source algorithms used to mitigate ray effects are summarized in Section 4. Section 5 contains some verification results and analysis. Concluding remarks are presented in Section 6.
2. Discretization Methods
The Boltzmann transport equation solved in ARES isThe physical meaning of each parameter is shown as follows: is angular neutron flux; is scalar neutron flux; is neutron velocity; is energy; is angle; is macroscopic total cross section; is macroscopic scattering cross section; is energy spectrum for prompt neutron; is external source term.
As we can see, the state is determined by the angular flux . The independent variables are (in centimeters), (in megaelectronvolts), and (in steradians). The following boundary condition is obeyed:which defines the incoming flux on all problem boundaries with outgoing normal .
Each independent variable needs to be discretized in deterministic solutions. The most commonly used angular discretization is the discrete ordinates method. This discretization produces a large, sparse, and linear system of equations in a six-dimensional phase space. We now briefly discuss each of these discretizations.
2.1. Angular Discretization
In method, we use a set number of discrete directions to discretize the continuous angular variable. With this approximation, (1) becomes
Quadrature sets are comprised of discrete directions and associated weight coefficients. For 3D geometry, ARES provides the level symmetric quadrature sets, the equal weight quadrature sets, and the even-odd moment quadrature sets, which are full-symmetric, as well as the Legendre-Chebyshev quadrature sets which are half-symmetric. And a biasing technique based on the Legendre-Chebyshev quadrature sets, called “the angular refinement technique for polar angles,” is developed.
For 3D transport calculation, the level symmetric quadrature sets (LS or ) developed by Lathrop and Carlson are still the most widely and commonly used sets [11]. level symmetric quadrature sets are illustrated by Figure 1. Nevertheless, the order of this type is limited to since weights become negative when we solve equations of moment conditions with the order is greater than .
From the point of view to satisfy more moment conditions for integration accuracy, the even-odd moment quadrature sets () have the ability to accurately integrate not only even moments but also odd moments of direction cosines in one octant [12]. On the other hand, the arrangement of directions in the sets may be more reasonable, considering that these sets have more degrees of freedom for direction cosines than LS and sets. Similar to LS sets, sets also have a limitation that their order cannot be greater than .
The half-symmetric Legendre-Chebyshev quadrature sets () are offered for users, whose set of direction cosines along axis is different from the sets along other two axes [13]. Legendre-Chebyshev quadrature sets are illustrated by Figure 2. Similar to LS sets, directions in sets are arranged along levels, equal to the roots of the Legendre polynomial or the one-dimensional Gauss-Legendre quadrature sets. The azimuthal angles or the direction cosines of and for each direction are defined from the roots of Chebyshev polynomial of first kind. From the Gauss-Legendre quadrature sets, the level weights for levels are determined and weights of directions on the same level are equal. The order of the sets can be easily increased without limitation of negative weights.
For problems that angular flux distribution or spatial distribution of scalar flux is highly peaked caused by highly anisotropic scattering or regional materials, results of transport calculation with traditional quadrature sets may be unsatisfying. To solve this problem, we designed a quadrature’s biasing technique called “the angular refinement technique for polar angles” based on the sets. We can refine arbitrary level, which users choose, to several levels to increase calculation accuracy in the local region. Especially, this technique can refine many levels simultaneously whether they are next to each other or not.
2.2. Spatial Differencing
ARES currently provides diamond difference with or without linear-zero flux fixup, theta weighted (TW), directional theta weighted (DTW), exponential directional weighted (EDW), and linear discontinuous finite element spatial differencing schemes.
The diamond difference method, which assumes a linear relationship between the directional flux at the cell center and cell boundaries, is simple and accurate for small mesh intervals. When the mesh interval is too large, the difference equations may yield negative fluxes. The TW, DTW, and EDW variations on the DD method were developed to eliminate the appearance of negative fluxes without significantly sacrificing computational cost or accuracy.
The balance equation can be obtained by integrating the discretized form of the transport equation over the mesh cell such as .where is the cell average flux and the entering and exiting angular fluxes are referred to as the “in” and “out” subscripts. is known from the previous source iteration and the entering angular fluxes are known from the boundary values. where ; ; .
The accuracy of the diamond difference scheme is second-order truncation. Considering that negative boundary angular flux is nonphysical, the negative flux set to zero fixup is commonly used. However, the fixup causes DZ to become nonlinear and depart from second-order accuracy [15]. The oscillation of DZ difference scheme is still apparent even when the mesh is refined. To guarantee a nonnegative exiting flux value with positive sources, the theta weighted (TW) scheme is developed.
The scheme uses the incoming fluxes to calculate weighting factors. The cell-centered and exiting fluxes vary smoothly between the step and diamond difference approximations. The weighting factors are calculated from the following system of equations:Using the weighting factors, the cell-centered and exiting fluxes are
The theta-weighting factors and are set to values between 0 and 1. By default, ARES uses the theta weighted model from the TORT code in which [16]. The weighted factors are bounded between the diamond difference and step approximations, .
The directional theta weighted scheme is an extension of TW scheme. The direction-based parameters are used to obtain angular flux weighting factor. To be consistent, the weights () are restricted to the range between 0.5 and 1. Because of the directional weighting of DTW, the over- and underestimated angular fluxes among different directions result in more accurate scalar fluxes [17]. However, it may be not a highly accurate scheme in all situations such as streaming problems. Therefore, the EDW scheme is developed.
The EDW scheme uses the DTW to predict a solution that is then corrected by an exponential fit and it should be more stable and accurate than DTW alone. We simply write down the equations that are solved in ARES. The inherently positive exponential auxiliary equations are given. The exponential coefficients can be obtained from the DTW solutions.where
This method is absolutely positive, stable, and directionally weighted and is significantly more accurate than the DTW scheme in streaming problems with relaxed cell intervals [18].
Discontinuous finite element differencing captures discontinuities in solution and material properties and has third-order accuracy for global quantities, is acceleratable and damped, and has the diffusion limit [19]. The DFEM spatial differencing remains accurate and stable even on coarse meshes [20]. We begin the discussion by assuming that the problem domain has been divided into unstructured tetrahedral volume elements. The material properties within each element are assumed to be constant. Carrying out the Galerkin finite element approximation [21] in which , equation for element becomeswherewhere is the inflow boundary, is the outflow boundary, and is the neighbored element of .
3. Solution Techniques and Parallel Algorithms
3.1. Numerical Solution Techniques
ARES uses traditional source iteration and Krylov methods to solve transport equation. To make the discussion of numerical solution techniques clear, discretized transport equation can be expressed in operator notation [22].where is the transport operator, is the operator that converts harmonic moments into discrete angles, is the discrete to moments operator that integrates discrete angles into angular flux moments using quadrature rules, is the scattering matrix, is the fission matrix, is the block matrix fission spectrum, and is the block matrix of nufission cross section.
The standard way to calculate eigenvalue is power iteration [23]. Power iteration proceeds as follows:
Within each power iteration, the method for solving (13) is equivalent to that for fixed-source calculations:where represents external fixed source for fixed-source problems and fission source for eigenvalue problems.
More specifically, with the groups defined over the energy range , (15) can be decomposed into a series of coupled single-group equations. Gauss-Seidel iteration is commonly used over energy and can be written as follows:
When using Gauss-Seidel iteration over energy, one must solve within-group equations over angle-space, and they have the following form:where involves all sources for group except for sources scattering from group .
ARES provides source iteration and Krylov for the within-group equations. Source iteration can be thought of as a two-part process:
Here, is the iteration index. It begins with an estimate flux moments for the scattering source on the right side of (18). By a transport sweep, new flux moments can be given through (18) and (19). These iterations are repeated until the flux moments converge.
However, as the problem becomes more scattering-dominated, source iteration will be increasingly inefficient. Classic diffusion synthetic acceleration scheme suffers from severe stability problems in three dimensions with large material discontinuities [24]. Therefore, except for source iteration, Krylov iterative method preconditioned with DSA is applied to solve within-group equations.
The desired form for Krylov iteration is given as follows:where
Krylov iteration schemes are particularly amenable for this quite large, fairly sparse matrix because only the action of operator on an iteration vector is required. Applying the action of operator on iteration vector requires some steps [25]. The transport operator is inverted by sweeping through the mesh in the discretized direction of particles travel. GMRES(m) is used as the Krylov solver in ARES code system.
3.2. Parallel Algorithms
The problems typically of interest in the nuclear engineering community are of large scale. As larger computer resources have made it possible, some discrete ordinates codes on massive parallel machines, such as PARTISN and Denovo, have been developed. Motivated by the required ability to calculate large scale problems on available computer resources, ARES is developed with capability of paralleling.
The efficiency of the discrete ordinates method is largely dependent on the efficiency of the transport sweep procedure. KBA algorithm [26], the best known parallel sweep algorithm devised by Baker and Koch, is used in ARES transport code system to obtain the ability of calculating large-scale problems. The KBA algorithm uses special columnar domain decomposition and a particular sweep ordering. Figure 3 depicts a typical cubic geometry, which has been divided into , , and meshes along -, -, and -axes, respectively.
For a given direction, the KBA algorithm orders the tasks as depicted in Figure 4. First, the processor that has been assigned task at the top front right corner solves the block. The solution of this block is depicted by removing from the mesh in Figure 4. Then, the newly computed partition boundary fluxes are sent to the neighbor processors along -, -axes direction. Along -axes direction, communication is not needed since all of these meshes belong to the same processor. Within each block, the sweep order is not important, as long as it satisfies the dependency constraints.
The discrete directions are pipelined so that the sweep along the next direction in the octant begins without waiting, when sweeping along the current direction is completed. The pipelining can be visualized as repeating the domains times, where is the number of discrete directions in an octant. Note that the factor of 2 results from the pipelining of octant within each quadrature.
We note that some processors remain idle at the beginning and ending of pipeline, which causes PCE is less than unity. Some more gains could be obtained by sweeping on all quadrants at once rather than sequentially [27]. This approach has not been included in current ARES code system.
4. Ray Effects Mitigation Methods
Ray effects, shown as unphysical oscillations in the scalar flux, are inherent problem in method. They are caused by the inability of any quadrature set to accurately integrate the angular flux. First collision source method [28] has been employed to mitigate ray effects in ARES.
The first collision source method decomposes the transport equation into (22) and (23) [28]. It means the total fluxes are decomposed into uncollided and collided fluxes.
The uncollided flux moments are calculated analytically by (24). The uncollided fluxes are applied to calculate first collision source with (25). The collided components are obtained with method in (23).
In ARES, arbitrary number of point sources are employed to approximate the source region. And the location of a point source can be arbitrary. RAY employs ray tracing method to accelerate the calculation of optical distances between all source points and all grid cells, and point source correction factor is introduced to improve the accuracy of calculation results [29]. RAY3D has been validated by a series of national benchmarks, such like Kobayashi benchmarks [30] and Azmy benchmarks [31]. RAY can effectively eliminate ray effects and obtain reasonable results.
5. Verification and Discussions
The ARES code system has been validated and verified by lots of analytical problems, international benchmark problems, and international benchmark experiments. All results of ARES have been compared with authoritative transport codes, such as TORT and MCNP.
The Kobayashi benchmark problem [30], proposed by OECD/NEA, can be used to verify the ability to calculate shielding problems with void region. The solutions of DONTRAN3D and DONTRAN3D with RAY3D agree with MCNP solutions that are within 3% and 11%, respectively [29]. ARES can effectively mitigate ray effects and obtain good accuracy for the void problems.
The Takeda benchmark problems [32], proposed by Takeda and Ikeda, are used to verify the accuracy of transport codes for critical calculation. The eigenvalues and region-averaged fluxes are calculated by DONTRAN3D, eigenvalue differences between ARES-calculated values, and reference values lie within 30 pcm in most cases [33]. ARES has a good performance in critical calculation for homogeneous core models.
In this chapter, the HBR-2 benchmark problem [34] and C5G7 benchmark problem [35] are calculated to further verify and validate the ability of ARES for engineering application.
5.1. HBR-2 Benchmark Problem
An accurate calculation of the neutron fluence at the reactor vessel is necessary to estimate the structural integrity over the designed lifetime and to support analyses for a potential plant life extension. The H. B. Robinson Unit 2 Pressure Vessel Benchmark (HBR-2 benchmark), based on experimental data from an operating PWR reactor, is the only openly available RPV benchmark through the SINBAD Database at the OECD/NEA Data Bank [14]. To verify the reliability and availability of ARES shielding calculation, the HBR-2 benchmark was modeled and we provided the final result of the average ratio of the calculated to measured specific activities () for the six dosimeters in the surveillance capsule during cycle 9.
HBR-2 benchmark is a 2300 MW (thermal power) pressurized water reactor designed by Westinghouse, as shown in Figure 5. The core consists of 157 fuel elements and is surrounded by the core baffle, core barrel, thermal shield, pressure vessel, and biological shield. The overall dimensions of the three-dimensional configuration is 443.31 × 443.31 × 425.936 cm, while each assembly is 21.504 × 21.504 cm, and each fuel assembly is composed of 15 × 15 array of fuel pins.
To describe the model accurately, the calculation of neutron source is significant. The power to neutron source conversion factor was calculated based on the contributions of 235U and 239Pu to the fission neutron source, and we took the average fission spectrum of 235U and 239Pu as the source energy spectrum [34].
The MUSE1.0 [36] cross-sectional library based on ENDF/B-VII was used for the ARES transport calculations. The approximation is introduced. corresponds to the order of the expansion in Legendre polynomials of the scattering cross-sectional matrix and represents the order of the flux angular discretization. Fully symmetrical quadrature sets were introduced.
With the purpose of comparing the calculated and measured specific activities, we take the reactor power changes during irradiation period and the affection of the closest fuel assemblies to the reaction rate at the dosimetry locations into account. According to the relative reactor data given by the benchmark, the approximate reaction rate can be written aswhere is reaction rates of the dosimeters in the surveillance capsule during th day, is reaction rate obtained from the transport calculation (ARES) in nominal core power, is normalized average relative power of the three closest fuel elements during th burn up step, is average relative power of the three closest fuel elements during the whole 9 cycles, is daily average reactor core power during th day (during th day is in the burnup step th), and is nominal core power (2300 MW).
For all the considered typical nuclide reaction, the reaction rates of the dosimeters in the surveillance capsule calculated by ARES, the reference values provided by DORT for the cycle average power distribution, and core power of 2300 MW are given in Table 1.
According to the reaction rates from Table 1, the specific activities were calculated and presented in Table 2. This table also lists the measured specific activities of the dosimeters from surveillance capsule at the end of cycle 9.
The ratios of the calculated and measured specific activities are listed in Table 3. The data in the bracket below the corresponding ratios are the values that take the correction discussed above into account.
The results indicate that the ARES transport calculation and the measured specific activities are in good agreement except reaction. However, the benchmark report pointed out that, due to the particularity of detectors, the deviation is generally large and further revision is needed. Overall, ARES has been proved that it can be applied in shielding calculation, and the results are relatively accurate and credible.
5.2. C5G7 Benchmark Problem
The C5G7 benchmark problem [35], proposed by OECD/NEA expert group, was modeled to verify the criticality calculation capabilities of ARES code system for reactor core problems without spatial homogenization. As shown in Figure 6, this benchmark consists of four UO2/MOX fuel assemblies surrounded by water moderator. Each fuel assembly is made up of a 17 × 17 square lattice of fuel pin cells. The pitch of each pin cell is 1.26 cm and every fuel pin, guide tube, control rod, and fission chamber have a radius of 0.54 cm. In our calculation, curved surface of the fuel rod is approximated by staircase. As shown in Figure 7, each of te pin cells is divided into meshes, which gives a total of meshes, together with moderator region. A detailed description of the fuel assembly composition is available [37]. Figure 7 illustrates that the volume of fuel and moderator is approximated by employing orthogonal structured meshes. Therefore, the macroscopic cross sections presented in [34] are corrected by adjusting nuclei densities to preserve fuel and moderator masses [37].where denotes actual volume and denotes approximate volume based on the orthogonal structure meshes.
This problem is executed using an EO16 quadrature set with a total of 288 angles. Three cases characterized by control rods position are considered: (1) UNRODDED, control rods stay in the moderator above fuel assemblies, (2) RODDED A, the control rods are inserted one-third of the inner UO2 fuel assembly, and (3) RODDED B, the control rods are inserted two-third of the inner UO2 fuel assembly and one-third of the two MOX fuel assemblies. The convergence criteria were set to for the flux and for the eigenvalue.
Table 4 gives the eigenvalue results and differences relative to the reference MCNP value. All of the three cases under consideration achieved excellent agreement with the reference solution within 20 pcm, calculated by
Figure 8 plots relative differences between ARES-calculated maximum pin power and reference MCNP values. Axially, the maximum relative differences lie in slice 3, in which region noticeable differences are presented with increased control rods insertion. This is likely due to the shortcomings of the method in treating large thermal flux gradients caused by inserted strong-absorptive control rods.
Assembly powers are calculated and compared with reference MCNP values in Table 5. The results indicate an overestimation of assembly power in the inner UO2 assembly. When moving toward the MOX assembly and outer UO2 assembly, it gives underestimated assembly powers. The differences all lie within one standard deviation uncertainty of MCNP calculations.
In summary, the results show that ARES can be used to perform precise transport calculation for complex three-dimensional geometries that include strongly absorbing materials, such as fuel assemblies without homogenization.
6. Conclusions
ARES code system is developed to solve the linearized Boltzmann transport equation for a wide variety of radiation transport applications and reactor physics analysis. ARES provides five spatial differencing schemes, uses the first collision source method to eliminate or mitigate the ray effects, and applies source iteration and Krylov iterative methods to solve the linear system of equations.
In this paper, ARES solutions to the HBR-2 benchmark and C5G7 benchmarks are in excellent agreement with reference results. ARES is undergoing continuous development with many new features planned for implementation, aiming to deal with advanced physical model and problems of extended range, such as more accurate solutions for deep penetration problems with void region.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by National Natural Science Foundation of China (11505059 and 11575061) and the Fundamental Research Funds for the Central Universities (2016MS59).