Corresponding author: Gennady Zherdev ( jerdev@ippe.ru ) Academic editor: Yury Korovin
© 2018 Gennady Zherdev, Tamara Kislitsyna, Mark Nikolayev.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Zherdev G, Kislitsyna T, Nikolayev M (2018) ROCOCO: A constants supply system for MonteCarlo reactor calculation. Nuclear Energy and Technology 4(2): 103109. https://doi.org/10.3897/nucet.4.30662

The ROCOCO system (
It should be noted that the system contains descriptions of detailed dependences of elastic scattering crosssections and angular distributions on all multiisotope elements; the relationship between the scattering angle and the energy loss in this case is determined with the use of the energydependent effective atomic weight.
The system’s programs are written in the FORTRAN language. The system is easily integrated in MonteCarlo codes.
Nuclear data, radiation field MonteCarlo calculation, combination of detailed and group descriptions of crosssection energy dependence
ROCOСO (Russian abbreviation for Calculation and Organization of COmbined COnstants) is a system designed to supply constants for MonteCarlo calculations of neutron fields and the photon fields they generate in nuclear reactors. It enables detailed description of the neutron crosssection energy dependences (as it is done in MCNP, a commonly known code, (
The ROCOCO system programs are written in the FORTRAN language. The system operability was tested using various compiler versions (GNU, LAHEY, INTEL, etc.). The system consists of three parts.
The first of these is the COMBINE package of programs designed to form a binary library of initial constants based on evaluated data files in the ENDF6 format (RUSFOND, a Russian library of evaluated neutron data files) and of group constants in the BNAB format (
The second part is used directly in the calculation. The procedures in this block read the computational assignment and an online library is formed on its basis for the required composition of materials and in the format convenient for simulation. In particular, inelastic scattering, thermalization, and fission spectra matrices are presented as accumulated probabilities of processes. It is also here that subgroup parameters of resonance isotopes are generated for the preset temperatures and compositions, while constants with detailed tracking of the neutron energy are recounted for the preset media temperatures. All this data is generated only for the isotopes present in the medium and is recorded into a dedicated exchange file that does not depend on the computation platform (what is meant is a binary or text file of a universal format read similarly by various compilers and in different operational systems, this enabling the user to prepare a library, e.g., on an office computer and to perform calculations on a supercomputer). The procedures in this block are united into an independently executed program module, COLIBRY (
The third part is intended to be used directly in MonteCarlo simulation programs and consists of several procedures. The key procedure, COLIBRY, is called only once in the process of a calculation. It refers to the exchange library prepared by the system’s second part, and accommodates them in respective general arrays which are accessible to all operating procedures.
The operating procedures are listed below. The first of them, GEFIS, has the parameters INUC, NGFIS, ig, E, and ves (here and hereafter, the parameters of the procedures are underlined for identification). Evidently, the fission spectrum is defined by which INUC nuclide fission gives birth to the neutron and which NGFIS energy group the neutron that has caused the fission belonged to. The procedure based on these parameters defines the group ig and the energy E of the fission neutron. The parameter ves is equal to unity by default. It is possible to increase artificially the yield of highenergy neutrons, in which case the parameter ves serves to compensate for the fission spectrum distortion. It is assumed in the first generation, when no INUC and NGFIS have yet been defined, that the neutron formed by the fission of the first of the listed neutrons (a priori fissionable) under the action of neutrons in the 50^{th} group (in which all fissile isotopes have fission crosssections other than zero). We shall note that there is an extended version of this procedure, called GEFISBETA, which, apart from the functions performed by the GEFIS procedure, forms the vesD array, eight components of which represent probabilities that the emitted fission neutron (in the group ig and with the energy Е) has been formed by the decay of the precursor for the respective delayed neutron group. The GEFISBETA procedure is intended for use in calculating effective fractions of delayed neutrons.
The SIGMA procedure has the parameters SIGMAC, SABS, SFNU, KERM, E, NG, and m. For a neutron with the energy Е, lying in the NG group in the physical zone m, it defines the full macroscopic crosssection SIGMAC, SABS is the macroscopic crosssection of absorption, the “fission neutron formation crosssection” SFNU = n×S_{f} , and KERM is the kerma value. The group’s energy and number are exchanged among the operating procedures automatically. To determine these values, the SIGMA procedure calculates initially the microscopic constants of all isotopes in the zone m for the energy of the considered neutron. If the neutron has formed in the considered zone, then the crosssections of all nuclides in this zone are determined for the energy of this neutron. If the neutron enters this area after crossing its boundaries, microconstants are determined only for the nuclides for which they have not been determined in the previous zones on the neutron flight path. Microconstants are entered in shared arrays to be further used in the COLLY procedure and, where required, in recording procedures (other than being part of the ROCOCO system). The SIGMA procedure also processes data for gamma quanta to which groups numbered 301 through 427 are attached (if the calculation assignment requires the gamma quanta paths generated by neutrons to be traced: the number of gamma groups is shown as equal to 127).
The COLLY procedure has the parameters m, NG, E, r, n, etha, NG1, E1, mulab, yi, ngam, and ygam. For the collision of a neutron taking place in the zone m, with the energy Е lying in the group NG, it defines the following type of collision: r = 1 is the elastic scattering; r = 2 is the inelastic scattering; and r = 3 is the absorption. Also n is determined as the number of the nuclide with which the collision took place; the average number of secondary fission neutrons formed in this absorption, E1 and NG1 are the energy and the number of the scattered neutron group; mulab is the scattering angle cosine in the laboratory coordinate system; yi is the set of scattered neutrons (exceeds unity if there is a contribution of the (n, xn) reactions to the inelastic scattering); ngam is the group number of the gamma quanta emitted in the inelastic scattering or capture; and ygam is the number of these quanta (301 ≤ ngam ≤ 427).
Therefore, it is enough for the connection of the ROCOCO system to the reactor MonteCarlo code that the COLIBRY subprogram (to the inlet of which only the exchange library file is supplied) is called prior to the solution of a particular problem. Then, the GEFIS procedure (defining the fission neutron energy and group) shall be called at the beginning of each selected path. The SIGMA procedure is called prior to tracing the track of the neutron (or the photon) lying within the zone in which the considered neutron or photon was born as the result of a process, or within the zone it entered when crossing the boundary. The COLLY procedure is called for the consideration of each collision.
Therefore, the procedures provide the calling program with all information required for tracing the neutron path and, where required, the path of the gamma quantum born in the neutron reaction. Also values are provided which are required to estimate the breeding ratio and the energy release in terms of range lengths, collisions and absorptions. Since, for each of the neutron track portion lying within the physical zone, the SIGMA procedure determines microconstants for all of the nuclides it contains, the program referring to ROCOCO may include recording procedures that estimate such functionals as the isotope reaction rates, disabled multigroup constants, etc.
The following is shown in the input file for the second block procedures the COLIBRY module includes:
MZONE – the number of physical zones (i.e., of zones differing in composition or temperature);
NIS – the full number of nuclides in at least one of these zones;
the number of the BNAB neutron groups MGRN ≤ 299 – in groups with numbers larger than 227 (i.e., with energies smaller than 4.64 eV), the thermal motion of neutrons is taken into account in calculations; a smaller number of groups (e.g., 298) needs to be specified for taking this into account;
the MGRG gamma group number equal either to zero or to 127 (as stated above, group numbers 301 through 427 are attached to gamma quanta);
NUCLIST – a list of NIS nuclides: foursymbol names of isotopes and twosymbol names of multiisotope elements;
ID – the array of NIS signs for the level of the crosssection structure description detail for each nuclide: ID = 2 with a detailed description of the crosssection structure; ID = 1 for the use of a 299group approximation with regard for the intragroup crosssection structure by subgroup method; ID = 0, a 299group approximation without regard for resonance selfshielding;
THERM – the set of the numbers of thermalization matrices for the BNAB system of constants (if this number is assumed to be equal to zero, then thermalization will be taken into account in a free gas approximation);
TZONE – the array of the physical zone temperatures (MZONE of values) (nuclides with ID = 1 or 2 contained in zones with different temperatures will be considered as different nuclides in the calculation);
RO – a twodimensional array of nuclear concentrations: concentrations of all nuclides should be shown for each zone.
It should be noted that the above keywords of the assignment parameters can be substituted by analogs adopted in the CONSYST constants supply system (Nikolayev M.N., VANT. Ser.: Yadernye konstanty. 1999, No. 2, p. 148).
In calculations, the crosssections of “detailed” nuclides outside the resonance region are described with the same level of detail as they are presented in evaluated data files. Detailed energy dependences of the crosssections, calculated for the temperature of 300 K using the RECONR and BROADER modules of the NJOY 2В code, are used in the region of resolved resonances. For zones with temperatures higher than the indoor temperature, detailed crosssections are recalculated using the DOPPLER procedure that implements the algorithm similar to that used in the BROADER module (
The crosssections of “subgroup” (ID = 1) and group (ID = 0) nuclides are described in a 299group BNAB approximation. And the resonance structure of crosssections for “subgroup” nuclides in the region of both unresolved and resolved resonances is described using subgroup parameters obtained based on tables of selfshielding factors interpolated to the temperature of the zone under consideration. In the region of unresolved resonances where two subgroups are normally enough to describe the structure of crosssections, subgroup parameters are determined using a simplified algorithm. In the region of resolved resonances where two subgroups are often not enough, the universal procedure contained in COLIBRY is used to obtain subgroup parameters based on selfshielding factors (
Energyangular distributions of neutrons in the thermalization region are described either in a multigroup approximation using the BNABRF thermalization matrices, or based on a free gas model using the algorithm described in (
The energy release in conditions of elastic or inelastic scattering or during entrainment or fission is calculated using 299group BNAB constants. The energy dependence of energy release within each group has the same level of detail as the description of crosssections. Where the calculation assignment requires the distribution of gamma quanta formed in neutron reactions (MGRG = 127) to be calculated, tables of local energy release are used, or, in the opposite case, tables of full energy release are used which also take into account the energy entrapped by gamma quanta. In both cases, the energy release data takes into account the contribution from the decay of neutron reaction products (with a halflife of up to three years).
Formation of gamma quanta is taken into account for the entrainment, inelastic scattering and fission reactions. 299×127group matrices of gamma quanta yields from the BNABRF are used. These matrices also take into account the yield of gamma quanta in the decay of products from respective reactions. The gamma quanta transport is calculated in a 127group approximation using macroscopic constants. The anisotropy of both coherent and incoherent scattering is taken into account by specifying probabilities of scattering for each of the three specified angles with respect to each intergroup transfer. The probabilities and the angles have been determined from the condition of five angular momentums being preserved.
We shall note the specific nature of the neutron interaction with nuclei of natural multiisotope elements. The format of the evaluated data file libraries enables representation of data for natural mixtures of multiisotope elements, but, in reality, this capability is used to describe the crosssections of only several materials rarely employed in engineering practice (e.g. such as platinum). ROCOCO enables a detailed description of the crosssection structure for all natural multiisotope elements. The crosssections of all isotopes are reduced to a single energy mesh which is thinned by dropping the energy nodes in which the crosssections of elements can be restored by linear interpolation from neighboring nodes with the preset accuracy (normally 0.1%). Energy dependences of the crosssections of elements are determined by adding the isotope crosssections to the weight of their concentrations. Angular distributions of neutrons, elastically scattered on isotopes and presented in the form of expansions in Legendre polynomials, are taken directly from the evaluated data files. The expansion coefficients for all isotopes are reduced to a single energy mesh and are averaged then with the weight of concentrations and the elastic scattering crosssections. At each point of the energy mesh, detailed dependence of the scattering probability on the angle of scattering in the inertia center system is restored after which 33 boundaries of 32 equiprobable intervals are determined from the scattering angle. In parallel with the crosssections, the energy dependence of the effective element mass, determined such that the average loss of energy during scattering is preserved, is calculated and stored. As in the event of individual elements, energyangular distributions of inelastic scattering are described in a multigroup BNABRF approximation. The use of the energy dependences of crosssections for natural elements reduces greatly the number of the nuclides the medium contains, which reduces the cost of the machine time for the calculation. The drift of the results when using this approximation is practically negligible.
To check the serviceability of the ROCOCO system, a program, kinfin, was written to calculate the criticality of infinite media. It enables estimation of breeding ratios for absorptions, k_{inf.0}, and for collisions, k_{inf.1}. No path lengths were estimated since, in an infinite homogeneous medium, a flying neutron will definitely experience a collision in this medium, and the estimates for collisions and path lengths will differ only in that the error will be somewhat higher in the latter case (the sample estimation error for the average path length will add to the estimation error for collisions). The multiplication of neutrons thanks to the (n, xn) reactions was taken into account by multiplying the neutron weight by the average multiplicity of the neutrons emitted as the result of inelastic scattering processes. The paths were calculated with the analog consideration of absorptions. The estimates obtained in the several first generations were not taken into account for all generations in the averaging of the results, as a steadystate spectrum of fission neutrons was established in these generations. After the selection among all generations, the results obtained for all generations (except a number of those dropped) were averaged, and the dispersions of each of the estimates (D_{0} and D_{1}) and the coefficient of correlation between them were computed (r). After that, based on the least square method, a combined estimate of the breeding ratio (k_{inf}) and its dispersion (D) were calculated:
comb = (1 – C)×k_{inf.0} + C×k_{inf.1},
where
The kinfin program was intended for verification calculations. The verification was done by comparing the calculation results with the data obtained using the MCNP5 code “tied” to the ROSFOND library. To make it possible to judge on the significance of the discrepancies in the results, media were chosen for the comparison, for which the results of the experimental k_{inf} determination were known (with some uncertainties). Another measure of the significance of the discrepancies was the difference between the results of the MCNP5 calculations using the ROSFOND and ENDF/B7 libraries. The MCNP5 calculations were performed in the same way as the kinfin calculations, with the analog consideration of absorptions. The results of the k_{inf} calculation for the four media will be presented below.
The first of these was uranium medium with an enrichment of 5.56%, SCHERZO5.56 (
kinfin: k_{inf} = 0.99750 ± 0.00041 (the counting time is 107 s),
MCNP(RF): k_{inf} = 0.99778 ± 0.00043 (the counting time is 174 s),
MCNP(B7): k_{inf} = 0.99288 ± 0.00043 (the counting time is 140 s).
Therefore, the difference between the kinfin and MCNP calculation results, amounting to 0.028%, is much less than the statistical error of this difference equal to ±0.56%. The counting time in the event of kinfin for this problem turned out to be one and a half time less than in the event of MCNP (obviously, thanks largely to the advantage of a multigroup description of the secondary neutron spectra from inelastic scattering and the (n, xn) reactions).
The results of the calculations based on the ROSFOND library agree with the experimentally estimated value of 1.0000 ± 0.0025 within the latter’s error limits.
The result of the calculation based on the ENDF/B7 library is 0.5% as less as the experimental value, which is twice as high as the experimental error.
A stainlesssteel medium with a small addition of highenrichment uranium was chosen for another test calculation and tested at the KOBRA test facility at IPPE (KBR9 assembly). The experiment conditions and the procedures for introducing computational corrections for the heterogeneity and the finite dimension of the investigated zone (having undergone an international expert appraisal) are described in detail in (NEA Nucl. Sci. Committee). The content of uranium with an enrichment of 90% was 90.2% by the number of nuclei. The rest of the components were iron, nickel, and chromium with small amounts of carbon, oxygen, silicon, titanium, and manganese impurities (Table
Composition of the KBR9 medium (numbers of the nuclei in the volume are presented, barn×cm).
C  3.1453E–4  Mn  8.2518E–4  Si  6.7256E–4  Ti  4.7341E–4  ^{235}U  3.7041E–4 
O  8.3126E–4  Fe  4.7217E–2  Ni  6.5656E–3  Cr  1.2788E–2  ^{238}U  4.0125E–5 
The experimental estimate of k_{inf} for this medium is 1.108 ± 0.009. As in the event of SCHERZO5.56, 1020 generations of 2000 paths each were selected among; the data for the first twenty paths were dropped. The results are as follows:
kinfin: k_{inf} =1.10807 ± 0.00051 (the counting time is 83 min),
kinfin (elem.) k_{inf} =1.10840 ± 0.00041 (the counting time is 64 min),
MCNP(RF): k_{inf} =1.10841 ± 0.00041 (the counting time is 103 min),
MCNP(B7): k_{inf} =1.11139 ± 0.00050 (the counting time is 102 min).
The result obtained based on kinfin, with the selection of paths by isotopes, coincided with the MCNP result obtained based also on the ROSFOND library. These results also coincide with the experimental value. Practically the same result was also obtained with the composition of structural materials defined by elements, which provides for a marked reduction in the counting time. The MCNP result, using the ENDF/B7 library, leads to k_{inf} being slightly overestimated.
Similar results were obtained also in calculations for the KOBRAinvestigated media with k_{inf} ~ 1 containing nickel (KBR7), molybdenum (KBR10), chromium (KBR15), and zirconium (KBR16) as the base material. The results of the calculations with the neutron energy traced in detail coincided with those obtained using MCNPROSFOND in the limits of statistical errors, but required a 20% smaller time. With the composition of the media having been described by elements, the time gain was 30 to 40% greater.
Two more test problems were considered. The first of them considered an infinite medium with the same composition as the unburnt core of the SVBR fast reactor. It contained 23 nuclides (including five fissile nuclides). The obtained results are as follows:
kinfin: k_{inf} = 1.2465 ± 0.0005 (the counting time is 20 min),
MCNP: k_{inf} = 1.2459 ± 0.0005 (the counting time is 21 min).
The counting time decreases to 17 min with the zone composition being described by elements.
The next problem considered the medium with the same composition as the burntup (to 11% of heavy atoms) SVBR reactor core containing 169 nuclides. The kinfin calculations described in details 52 nuclides the concentration of which exceeded 0.5%. The rest were described in a group approximation. The results are as follows:
kinfin: k_{inf} = 1.1436 ± 0.0003 (the counting time is 52 min),
MCNP: k_{inf} = 1.1437 ± 0.0003 (the counting time is 209 min).
The obtained data shows graphically the advantage of a group description of nuclides contained in miniscule concentrations. This is only natural since one has to search for the crosssections matching the energy resulting from the selection process for all nuclides independent of their concentration.
Finally, aqueous solutions of highenrichment uranylnitrate and plutonium were calculated. The characteristics of the solutions and the k_{eff} calculation results are presented in Table
Results of calculating the criticality of aqueous plutonium and highenrichment uranium solutions using different codes and libraries of constants
Solution  MMKRF  MCNP5 (ROSFOND)  MCNP5 (ENDF/BVII)  Difference from MCNP5 (ROSFOND)  
MMKRF  ENDF/BVII  
Pu+H_{2}O  0.99388(8)  0.99369(10)  0.99404(10)  0.02%  0.04% 
U+H_{2}O  0.90857(8)  0.90880(10)  0.90889(10)  –0.03%  0.01% 
No major differences in neutron spectra have been detected (it is only with energies of < 10^{–3} eV that the MNCP5 neutron fluxes turned out to be 2 to 5% smaller than those calculated based on kinfin).
Testing has shown that ROCOCO makes it possible to obtain the same results as the MCNP code with smaller time expenditures thanks to the group description of inelastic scattering indicatrices and with the possibility of a group or subgroup description of the crosssections of nuclides with a low concentration and a detailed description of the crosssections of elements (for which there are no evaluated data files).
The connection by V.B. Polevoy of ROCOCO to the MMKFK code (
It makes sense to use the ROCOCO system in precision codes for engineering calculations of power reactors in which increased speed of response will make it possible to solve problems more complex than criticality problems (solution of space and time problems, estimation of local functionals, etc.), to improve the quality of calculations, and to expand the application scope (
The authors express their profound gratitude to V.N. Koshcheyev, G.B. Lomakin, and V.V. Korobeynikov for their assistance with the study; to V.B. Polevoy who demonstrated the feasibility of connecting the ROCOCO system to engineering codes, and to Ye.V. Zhemchugov who developed and provided the SUBGGROUPS program, one of the code’s key subprograms.