79urn:lsid:arphahub.com:pub:D015427E-53F6-553C-9279-2E51CD684756Nuclear Energy and TechnologyNUCET2452-3038National Research Nuclear University MEPhI (Moscow Engineering Physics Institute)10.3897/nucet.9.102507102507Research ArticleReactor physicsInfluence of the spatial grid type on the result of calculating the neutron fields in the nuclear power plant shieldingNikolaevaOlga V.nika@kiam.ru1GaifulinSergey A.1BassLeonid P.1DmitrievDenis V.2NikolaevAlexandr A.3Keldysh Institute of Applied Mathematics, Russian Academy of Sciences, 4, Miusskaya Sq., 125047 Moscow, RussiaKeldysh Institute of Applied Mathematics, Russian Academy of SciencesMoscowRussiaIPPE JSC, 1 Bondarenko Sq., 249033 Obninsk, Kaluga Reg., RussiaIPPE JSCObninskRussiaOKB Gidropress JSC, 21 Ordzhonikidze Str., 142103 Podolsk, RussiaOKB Gidropress JSCPodolskRussia
Corresponding author: Olga V. Nikolaeva (nika@kiam.ru)
Academic editor: Yury Korovin
2023200620239299105D96A019C-B5F3-50A3-BAE8-E94A88896A6B1406202201032023This is an open access article distributed under the terms of the CC0 Public Domain Dedication.
The paper considers the influence of the spatial grid type on the result of solving the equation of neutron transport in the nuclear power plant (NPP) shielding.
Neutron fields were calculated in a realistic model of a liquid metal cooled fast neutron tank reactor with an integral equipment layout. Structured cubic and unstructured hexahedral grids (pmsnsys and FRIGATE codes) and unstructured tetrahedral and prismatic grids (RADUGA T code) are used. Limiting values of the group fluxes averaged over the material zones for refined grids have been obtained.
It has been shown that the calculation results depend on the type of approximation for the curvilinear inner boundaries between the material zones rather than on the grid cell type (cube, hexahedron, tetrahedron, prism). Using “toothed” approximations for curvilinear boundaries leads to an increase in the area of the boundaries, as well as to the neutron flux refraction condition arising on them. These effects lead to an upward bias in the transport equation solution, and for all energy groups.
Conclusion. When solving an equation of neutron transport in the NPP shielding by a grid technique, it is necessary to use grids other than leading to “toothed” approximations of the inner boundaries. Tetrahedral or prismatic grids, or grids of arbitrary hexahedrons can be recommended, as well as composite grids in which cubic cells are located inside the material zone, and hexahedron cells are located near the zone boundary.
transport equationgridscurvilinear boundariesCitation
Nikolaeva OV, Gaifulin SA, Bass LP, Dmitriev DV, Nikolaev AA (2023) Influence of the spatial grid type on the result of calculating the neutron fields in the nuclear power plant shielding. Nuclear Energy and Technology 9(2): 99–105. https://doi.org/10.3897/nucet.9.102507
Introduction
Important in solving the problem of neutron transport in the nuclear power plant (NPP) cores and shielding is spatial grid selection. In particular, the grid cell boundaries are expected to approximate, as accurately as possible, the boundaries of the materials. The earliest software packages proposed that calculations be performed in a Cartesian xyz geometry based on parallelepiped cells. To take into account the curvilinear boundaries, packages were developed to allow calculations in curvilinear rzϕ and rθϕ geometries. However, curvilinear geometries do not allow one to describe at all times in an adequate manner all boundaries between the materials. As results, calculations have to be run in Cartesian xyz coordinates (Alcouffe 1990; Azmy 1996). Further, this group of codes evolved towards vectored calculations (Dahl 2006; Humbert 2006; Corau and Sjoden 2011) and improved grid schemes and iterative methods for solving systems of finite-difference equations (CNCSN. One, Two- and Three-Dimensional Coupled Neutral and Charged Particle Discrete Ordinates Parallel Multi-Threaded Code System 2009; Moryakov 2017; Royston et al. 2018).
Using a grid of parallelepiped cells in xyz geometry leads to replacement of curvilinear boundaries as well as boundaries that are not parallel to the coordinate axes by “toothed” surfaces (we shall refer to a surface composed of the grid cell faces as a “toothed” surface if the angle between adjoining faces is potentially more than 180^{о} (see Fig. 1). It is proposed that calculations are performed using a very fine grid (with a cell size of ~ 2 cm (Kovalishin et al. 2018) when the “toothed” boundary is geometrically close to the real one.
Directions of the escaping and retuning neutrons for two types of grid boundaries where d is the segment length of the broken line that approximates the circumference in the cylinder cross-section; β is the angle between the broken line segments; 1 is the direction of the escaping neutron; 2 is the direction of the neutron that has escaped, returned and escaped again.
https://binary.pensoft.net/fig/866039
Better quality of the grid approximation for the boundaries of the materials requires hexahedron cells (PC Software. Calculation of Neutronic Characteristics of Nuclear Reactors. “PMSNSYS”. 8624607.00622. 2011). The further step is to use unstructured grids of tetrahedrons and prisms (Vassiliev et al. 2008; Aristova and Astafurov 2017; Chen et al. 2017; Kim and Lee 2017; Belousov et al. 2018; Hoagland et al. 2021; Nikolaeva et al. 2021), as well as grids of arbitrary hexahedrons (Nikolaev 2016). In this case, it becomes possible to accurately set all straight-line boundaries; curvilinear boundaries are approximated by “faceted” surfaces (we shall refer to the surface composed of the grid cell faces as a “faceted” surface if the angle between the adjoining faces does not exceed 180°, see Fig. 1).
Another approach is to use structured grids with parallelepiped cells with defining additional “composite” substances in the cells where the boundary pass through. The interaction cross-sections for each additional substance are calculated as a linear combination of the cross-sections within the respective cell of materials; the linear combination coefficients are assumed to be the mass fractions of each material in the cell (the so-called Volume Fraction method). A foreign grid generator (Orsi 2006) and a similar Russian grid generator (Gurevich et al. 2007) are known which offer such approach. Many codes include the capability for calculations with additional substances (Dahl 2006; Humbert 2006; CNCSN.One, Two- and Three-Dimensional Coupled Neutral and Charged Particle Discrete Ordinates Parallel Multi-Threaded Code System 2009; Royston et al. 2018). The mass fractions of materials are however just approximately equal to the required linear combination coefficients. It is therefore critical to generate a grid that defines the boundaries of materials as accurately as possible.
Models of the nuclear reactor cores and shielding normally include multiple cylinders. When structured grids of parallelepipeds are used, the cylinders are substituted for prisms with “toothed” side surfaces. When unstructured tetrahedral or prismatic grids or grids of arbitrary hexahedrons are used, the cylinders are substituted for prisms with “faceted” surfaces. By adjusting the position of the grid nodes, one can achieve the coincidence of the “grid” and the initial cylinder volumes. However, the area of the “toothed” and “faceted” prism surface is varied and differs from the initial cylinder surface area. Besides, when a “toothed” grid boundary is used, such effect takes place as the return of the escaping neutron back into the cylinder (see Fig. 1). It is essential to assess the effects of the grid type on the result of the grid simulation for the reactor shielding neutron fields.
Description of the grid model
A model of a lead-bismuth cooled fast neutron tank reactor with an integral equipment layout was developed for the study. The model comprises the core, the coolant (lead-bismuth), internals (coolant, steel), the evaporator (coolant, steel, water), and the water shield tank, and the whole of the structure is accommodated in a concrete container (Fig. 2). The model height is 731 cm, and the model radius is 400 cm. Only a fourth part of the region is considered with defining the reflection condition on surfaces x = 0 and y = 0 (see Fig. 2). A stationary source is defined in the core.
Model of a fast neutron tank reactor (planes z = 453 and y = 0): 1 – core; 2 – steel; 3 – coolant; 4 – in-tank shielding; 5 – evaporator; 6 – gas gap; 7 – water shield tank; 8 – concrete; 9 – air; 10 – steel and coolant; 11 – boron carbide; 12 – reflection condition; 13 – flux density output region (see Fig. 7).
https://binary.pensoft.net/fig/866040
Each material in the model is matched by one or more cylinders or cylindrical layers with axes that are parallel to the problem’s z axis (see Fig. 2). The following methods are used to generate spatial grids.
A uniform grid with step
r is introduced for each variable
x,
y, and
z. A regular grid of cube cells with the edge length
r is produced as the result. We shall denote such grid type by
C.
A grid of arbitrary hexahedrons with the upper and lower faces that are parallel to plane
z = 0 is generated using the REBEL‑III code (Dedul and Nikolaev 2014). Some of the grid cells are connected through their surfaces without joining at the vertices. The hexahedrons take the form of cubes far from the boundaries. We shall denote such grid type by
H.
Each circle in plane
xy is substituted for a regular polygon with the edge length
d; as a result, each cylinder is substituted for a faceted prism. The Salome code (https://www.salome‑platform.org) is used in this region to generate an unstructured tetrahedral grid with the mean cell radius
r (the cell radius is the radius of a ball of the same volume). We shall denote such grid by
Tf.
Each circle in plane
xy is substituted for a set of squares with the side length
d; as a result, each cylinder is substituted for a toothed prism. Further, an unstructured tetrahedral grid with the mean cell radius
r is generated in the obtained region using the Salome code. We shall denote such grid by
Ts.
An unstructured triangular grid with the mean cell edge length
r is generated in plane
xy using the Salome code. A grid with step
r is introduced axially. As a result, an unstructured grid of triangular prims is produced. We shall denote such grid type by P.
Each grid (Table 1) is characterized by the mean cell radius (step) r and the segment length d of the broken line that approximates the circumference in the cylinder cross-section (Fig. 3). Parameter d is defined a priori for grid types Tf and Ts. For grid type P, quantity d is defined as the mean length of the grid edges that form the evaporator boundary in plane xy. We assume that d = r for the rest of the grid types.
Grid approximations of the evaporator cross-sections: 1 – Tf, d = 25; 2 – Tf, d = 12.5; 3 – Tf, d = 6.25; 4 – P, d = 3; 5 – Ts, d = 25; 6 – Ts, d = 12.5; 7 – H, d = 1.6; 8 – C, d = 1.
https://binary.pensoft.net/fig/866041
Spatial grids used
Grid type
Cell type
Type of grid boundary approximating cylindrical surface
C
Cube
Toothed
H
Hexahedron
Faceted
Tf
Tetrahedron
Faceted
Ts
Tetrahedron
Toothed
P
Triangular prism
Faceted
Fig. 3 presents grid approximations for the evaporator cross-sections. Grids of types Tf, P and H form “faceted” approximations for the sides surfaces of the cylinders. Type Ts and С grids form “toothed” approximations. The grid nodes are selected such that the volumes of the “grid” cylinders are equal to the volumes of the original cylinders. No VolumeFraction method (Orsi 2006; Gurevich et al. 2007) is used here to keep the mass balance.
For simulations, we shall use uniform Carlson grid S12 (Carlson 1976) and ENDF/B‑VI.8 constants (30 groups of neutrons, energy interval (10^{-4} eV, 17 MeV)) (Voronkov et al. 2009).
Calculations on cubic grids of type C are run by the PMSNSYS code (PC Software. Calculation of Neutronic Characteristics of Nuclear Reactors. «PMSNSYS». 8624607.00622 2011) by means of the diamond DD0 scheme with correction of the solution for positivity. The FRIGATE code (Nikolaev 2016) and the DDL scheme with the correction of the solution for positivity is used for the hexahedral grid calculations. The RADUGA T code (Nikolaeva et al. 2021) with the finite element method is used for calculations based on unstructured grids of types Tf, Ts and P. The PMSNSYS calculations were performed on the SuperEVM computer at JSC OKB Gidropress, the FRIGATE calculations were performed on a multicore personal computer, and the RADUGA T calculations were performed on the K60 hybrid supercomputer installed at the Shared Supercomputer Center at the Keldysh Institute of Applied Mathematics of the Russian Academy of Sciences. Densities of the group neutron fluxes are found, and their average values in each of the material zones are computed.
Calculation results
Initially, different types of grids with the ever-decreasing cell radius r (the step for structured grids) are used for the calculations. Grids of types Tf and P are used with the radius values 4.7, 3.5, 2.4, 1.7, 1.5, and 1.2 cm. Type C grids are used with steps of 2 and 1 cm. Values Ψ obtained in the calculations for the evaporator-average (see Fig. 2) neutron flux density in group 1 are shown in Fig. 4. The convergence of values Ψ at grid refining is achieved with the cell radius (step) value being r ≈ 1.5 cm. Value Ψ obtained in grid refining does not depend on the type of the grid scheme used.
Average neutron flux densities in group 1 in the evaporator obtained in the process of refining grids of different types: 1 – C grids; 2 – P grids; 3 – Tf grids with d = 6.25 cm; 4 – Tf grids with d = 12.5 cm; 5 – Tf grids with d = 25 cm.
https://binary.pensoft.net/fig/866042
For each grid type however, when the cell radius (step) r decreases, values Ψ converge to different values (see Fig. 4). Limiting values Ψ for type Tf grids depend on the segment length d of the broken line that approximates the circumference in the cylinder cross-section (see Fig. 3). When the segment length d decreases, the limiting values obtained on type Tf grids converge to the limiting value obtained on type P grids. This value, however, does not coincide with the limiting value obtained on type C grids.
Fig. 5 presents limiting values Ψ obtained with different grids with the cell (step) radius being r ≈ 1.5 cm. It can be seen that the limiting value depends not on the cell type but on how the cylindrical boundaries of the materials are approximated (by “toothed” or “faceted” surfaces).
Average flux densities for neutron group 1 in the evaporator obtained based on grids with different cell types: 1 – “toothed” approximations of boundaries (C and Ts grids); 2 – “faceted” approximations of boundaries (Tf, P and H meshes). Letter symbols for cell types: с – cubes; t – tetrahedrons; p – triangular prisms; h – hexahedrons.
https://binary.pensoft.net/fig/866043
Let us find the limiting average values of the flux densities 〈Ψ〉 for different material zones in all energy groups. We shall perform the calculations on type C grids (“toothed” boundaries, r = 2 cm, 12 million cells) and type P grids (“faceted” boundaries, r = 1.7 cm, 4 million cells). We shall note that cells are much fewer when changing structured C grids for unstructured P grids. The calculations take quite a long time (for P grids – 76 hours, using 672 cores of the K60 supercomputer at the Keldysh Applied Mathematics Institute).
Fig. 6 presents ratios of values 〈Ψ〉 obtained on a P grid to the values obtained based on a C grid. The ratios are presented only for the material zones where they differ from unity by more than 0.01. It can be seen that the influence of the grid type is preserved in all energy groups though not across the region. Besides, it can be seen that the flux densities obtained based on a grid with “toothed” approximations are higher than the densities obtained on a grid with “faceted” approximations.
Ratio of multi-group neutron flux densities obtained based on a grid with “faceted” boundaries to the flux densities obtained based on a grid with “toothed” boundaries: 1 – evaporator; 2 – water shield tank; 3 – concrete.
https://binary.pensoft.net/fig/866044
This can be explained with the help of the diagram in Fig. 1. If the interaction cross-sections outside the cylinder are much smaller than in the cylinder, then the neutron that has escaped the cylinder through one boundary tooth may return into the cylinder through the neighboring tooth. In fact, the “toothed” boundary imposes the flux density refraction condition.
Fig. 7 presents the flux densities in group 1 on a part of the concrete container’s lower boundary (z = 0); the selected part of the plane is separated from the outer boundary by an air layer (see Fig. 2). The flux density obtained with “toothed” approximations of the boundaries is higher near the outer boundary than the flux density obtained with “faceted” approximations (see Fig. 7). As a result, the flux density ratio (see Fig. 6) is less than unity.
Flux densities in group 1 on the lower boundary (z = 0) of the concrete cask (9 – air; 8 – concrete): a.C grid; b.Tf grid (logarithmic scale).
https://binary.pensoft.net/fig/866045
The flux density upward bias is not uniform across the region. The water shield tank is situated against the concrete container and the returning neutrons rapidly reach the tank. The evaporator is situated deep within the region but its upper part is in the air space (see Fig. 2). The neutrons that have passed through the concrete container enter the same part of the evaporator at once. The upper and lower boundaries of the container are parallel to plane z = 0. “Excessive” neutrons enter the container only through the “toothed” boundary that approximates surface r = 400 cm. Therefore, a major upward bias for the flux density near boundary r = 400 cm (see Fig. 7) leads to a minor upward bias in the container-average value.
On the other hand, the surface with a “toothed” boundary has a larger area than that with a “faceted” boundary. The lower part of the evaporator is situated near the source. More neutrons enter the evaporator through the “toothed” boundary than through the “faceted” boundary (Fig. 8). This is one more factor of the flux density upward bias in calculations with “toothed” approximations of the boundaries.
Flux densities in group 1 on the lower boundary of the evaporator a. C grid; b.Tf grid (logarithmic scale).
https://binary.pensoft.net/fig/866046Conclusion
A question has been considered concerning the influence of the spatial grid type on the result of solving a multi-group equation of neutron transport in the NPP shielding using a grid technique. The study was undertaken based on a realistic model of a heavy liquid metal cooled fast neutron tank reactor with an integral equipment layout. Structured grids with cubic cells and unstructured grids with hexahedron, tetrahedron and triangular prism cells were used.
The use of cubic cells gives birth to “toothed” approximations of the curvilinear inner boundaries (between the materials). On the other hand, “toothed” approximations are generated for the boundaries as such, and tetrahedral grids are generated then in the obtained region. The generation of tetrahedral, prismatic and hexahedral grids in the initial region leads to “faceted” approximations of the curvilinear boundaries. The grid nodes are adjusted such that the total volume of the cells that match each material is equal to the volume of this material in the initial model.
The values of the group flux densities, average for the material zones, are found. The values of the obtained quantities turn out to depend not on the grid cell type (cube, tetrahedron, prism, hexahedron) but on how this grid’s cell faces approximate the curvilinear boundaries (“toothed” or “faceted”).
“Faceted” approximations of the boundaries impose the refraction conditions on them; these conditions lead to the flux density redistribution among the zones. The refraction conditions on the boundary between the concrete container and air (in fact, on the region’s outer boundary) lead to a part of the neurons that have escaped from the region returning to it.
Besides, the areas of the “toothed” boundaries are larger than those of the “faceted” boundaries. “Toothed” approximations increase the surface for the neutrons that have escaped from the core to enter the adjoining zones.
As a result, there is an upward bias observed for the flux densities in all energy groups. This upward bias is maximum on the periphery (reaching potentially 10%). And this quantity has been found for the integral flux density values (average for the material zones). The upward bias for local values is expected to turn out greater.
Therefore, when solving a transport equation using a grid technique, it is important to use such grids that do not lead to “toothed” approximations of the curvilinear boundaries. Tetrahedral and prismatic grids can also be recommended. Composite grids can be used in which the cube cells are inside the material zone, and the hexahedron cells are near the zone boundary (as the type H grid used in the study).
ReferencesAlcouffeRE (1990) Three dimensional transport benchmark exercise using THREEDANT. Proc. Of meeting on three-dimensional neutron transport benchmarks. Los Alamos, New Mexico. https://inis.iaea.org/collection/NCLCollectionStore/_Public/22/031/22031958.pdf?r=1 [accessed Jun. 10, 2022]AristovaENAstafurovGO (2017) The second order short characteristics method for the solution of the transport equation on tetrahedral grid.AzmyYY (1996) The three-dimensional, discrete ordinates neutral particle transport code TORT: an overview. Proc. of OECD/NEA meeting on 3D deterministic radiation transport computer programs, feature, applications and perspectives. Paris, France.BelousovVIGrushinNASychugovaEPSeleznevEF (2018) Some results of verification of code ODETTA for shielding calculations. VANT. Ser.CarlsonBG (1976) A method of characteristics and other improvements in solutions methods for the transport equations.ChenYZhangBZhangLZhengJZhengYLiuC (2017) ARES: A parallel discrete ordinates transport code for radiation shielding applications and reactor physics analysis. Hindawi Science and Technology of Nuclear Installations, Article ID 2596727. https://doi.org/10.1155/2017/2596727CNCSN One, Two- and Three-Dimensional Coupled Neutral and Charged Particle Discrete Ordinates Parallel Multi-Threaded Code System (2009) RSICC code package CCC-726.CorauTSjodenG (2011) 3D neutron transport and HPC: A PWR full core calculation using PENTRAN SN code and IBM BLUEGENE/P computers.DahlJA (2006) PARTISN results for the OECD/NEA 3-D extension C5G7MOX Benchmark.DedulAVNikolaevAA (2014) REBEL – software of pre and post processing calculations of the neutronic performances of reactor on high-velocity neutrons with lead bismuthic heat carrier.GurevichMIRusskovAAVoloschenkoAM (2007) ConDat 1.0 – code for converting by the tracing algorithm the combinatorial geometry presentation to the bit-mapped one. Users Guide. Preprint No 12. Keldysh Institute of Applied Mathematics, Russian Academy of Sciences Publ. https://keldysh.ru/papers/2007/prep12/prep2007_12.html [accessed Jun. 10, 2022] [in Russian]HoaglandDSYessayanRAAzmyYY (2021) Solution of the neutron transport equation on unstructured grids using the parallel block jacobi-integral transport matrix method via the novel green’s function ITMM construction algorithm on massively parallel computer systems.HumbertP (2006) Results for the C5G7 3-D extension benchmark using the discrete ordinated code PANDA.KimJWLeeYO (2017) A deep penetration problem calculation using AETIUS: an easy modeling discrete ordinates transport code using unstructured tetrahedral mesh, shared memory parallel. EPJ Web of Conferences 153: 06025. https://doi.org/10.1051/epjconf/201715306025KovalishinAAMoryakovAVRaskachKF (2018) Neutronics calculation of fast reactor using modern computing systems.MoryakovAV (2017) LUCKY_TD code for solving the time-dependent transport equation with the use of parallel computations.NikolaevAA (2016) Generalization of two-dimensional DDL schemes of the GQ method for three-dimensional arbitrary hexahedral spatial mesh.NikolaevaOVGaifulinSABassLP (2021) Code RADUGA T for simulating neutrons fluxes in nuclear power stations. Vestnik YuUrGU. Ser.OrsiR (2006) A general method of conserving mass in complex geometry simulations on mesh grids and its implementation in BOT3P5.0.PC Software. Calculation of Neutronic Characteristics of Nuclear Reactors. “PMSNSYS”. 8624607.00622 (2011) OKB “GIDROPRESS”. [in Russian]RoystonKEJohnsonSREvansTMMosherSWNaishJKosandB (2018) JET contributors. Application of the denovo discrete ordinates radiation transport code to large-scale fusion neutronics.VassilievONWareingTADavisIMMcGheeJBarnettDHortonJLGiffordKFaillaGTittUMourtadaF (2008) Feasibility of a multigroup deterministic solution method for 3D radiotherapy dose calculations.VoronkovAVSinitsaVVDedulAVKalchenkoVV (2009) Libraries of multigroup constants for Code REACTOR-GP. VANT. Ser.
Russian text published: Izvestiya vuzov. Yadernaya Energetika (ISSN 0204-3327), 2022, n. 3, pp. 146–157.