Computer models of gas-dynamic processes of radionuclides behavior in thermal emission converter of nuclear power plant

The paper is devoted to the application of numerical solution of the Boltzmann equation method for modeling the behavior of radionuclides (Xe and Kr) in the inter-electrode gap and the vacuum-cesium channel of a thermionic reactor-converter, which is used in a space nuclear power generator. In this work, we developed modeling techniques and implemented them in a software package, that is using the algorithm of simultaneous calculations for the left side of the Boltzmann equation. We performed series of calculations for different initial flux of radionuclides and evaluated the dependence of the radionulcides output flux on the input flux for two areas of the reactor.


Introduction
In the design development of space nuclear power plants, thermal emission reactor-converters based on thermal and intermediate neutrons were developed [1]. One of the elements of such a reactor is an electricity generating channel (EGC) containing a vented fuel element with nuclear fuel based on uranium dioxide. The operability of a multi-element EGC with oxide nuclear fuel with communicating cavities of a fuel rod and the interelectrode gap (IEG), and, in particular, a change in the output electrical characteristics is influenced by numerous factors. One of them is the presence in the IEG of fission products formed in fuel during reactor operation [2]. The fission metal products and their oxides enter the IEG diffusing through the trap and the exhaust path, condensing on the walls, however, the gaseous fission products, for example, Kr and Xe do not condense on the trap and the exhaust path, entering the IEG EGC and then the vacuum-cesium system [3]. In the process of diffusion of radionuclides through a cesium vacuum system, their beta decay occurs, resulting from Ba and Sr have a significant effect on the output electric power of EGCs due to adsorption on the walls of the collector [4]. To justify and refine the reactor circuit, tests are carried out in a loop channel [5], using a γ-spectrometric system, placed on a vacuum pipe of the pumping line, radionuclides are recorded and the volumetric activity of the pumped gases is measured. Promising carbon-based materials are described in the papers [6][7][8][9][10], relaxation processes in such materials are of particular interest [11][12][13][14][15][16][17].
A problem in the design of this installation is to determine the pressures and gas flows in the IEG and in the cesium working fluid steam generator (FSG), since estimates of the release of  radioactive gases from the fuel are time-consuming and do not allow determining the initial gas  flows from the fuel through an exhaust device with acceptable accuracy. To analyze the gaseous mixture in the EGC cavity, a numerical method is used to solve the Boltzmann equation [18], which describes the behavior of discharged gas mixtures. The Boltzmann numerical equation is local in space, which allows efficient calculations on multicore systems, performance losses occur only due to the exchange of information between neighboring cells distributed in different processes, space partitioning algorithms can reduce the number of neighboring cells and achieve a close to linear increase in productivity.

Device description
A multi-element EGC is a complex engineering structure that contains series-connected EGCs with an internal arrangement of fuel rods enclosed in an airtight envelope, switching adapters connecting emitters to collectors of neighboring power generating elements. Channels connecting the cavities of the IEG and the fuel element are made inside the switching adapters. The space of the IEG is filled with gaseous cesium -a working fluid, krypton and xenon enter the gap through the channels in the emitter; during diffusion, they decay along beta decay chains: Ba The estimated pressures Cs in the cavity of the IEG are 150 Pa, and Kr and Xe are in the cavity of the IEG under a pressure of no higher than 1.5 Pa. The characteristic dimensions of the IEG vary from 0.4 mm to 520 mm, this gives a spread of Knudsen number in the range from 0.004 to 100. This range allows us to consider the movement of gases in the gap as the movement of discharged gases and allows you to simulate their behavior by numerically solving the Boltzmann equation.
FSG consists of an evaporator, a pipeline and a condenser. The evaporator is two coaxial cylinders in the gap between which there is liquid cesium at a given temperature. The inner hollow cylinder has a hole system and a top cover. When cesium vapor enters the vacuum system, it condenses on the cooled walls of the pipeline (condenser) and flows into the evaporator, and non-condensable gases are pumped out by high vacuum pumps. Estimated pressures in the cavity of the FCG are 150 Pa for cesium gas, 10 −6 -10 −7 Pa for Xe and Kr, respectively. Temperature on the wall vary in the range from 320 K to 600 K, the characteristic dimensions of the installation vary in the range of 1-250 mm, which gives estimated Knudsen numbers for the gas mixture in the range from 0.0004 to 100. For such Knudsen numbers, the gas mixture can be considered sufficient discharged for modeling by the method of numerical solution of the Boltzmann equation.

Mathematical basis
To analyze the behavior of gases, the three-dimensional Boltzmann equation is used, and when calculating the collision integral, the approximation of the potential of hard spheres is used: where z, b, ε is a cylindrical coordinate system, the z axis of which passes through a molecule with a momentum p and is parallel to it, b is the impact distance (radius), and ε is the azimuthal angle, g is the relative velocity of approach, p, p 1 are the momenta of the molecules before the collision, p , p 1 are the momenta acquired by the molecules after the collision. To solve the system of Boltzmann equations, a monotonic difference scheme is used using the method of splitting according to physical processes [19]: at the first stage, the left side of the Boltzmann equation -the transport equation is solved, at the second stage the right side of the Boltzmann equation -the relaxation equation is solved, the processes of intermolecular interaction in gases: To solve the transport equation, a first-order accuracy scheme is used, the distribution function is specified in each spatial cell, the numerical form of the transport equation is: where τ is the time step, V i is the volume of the i-th spatial cell, k is the index of the face of the i-th cell, S ik is the face area between the i-th and k-th spatial cells, n ik is the normal to the face between the i-th kth cells. For each face, a boundary distribution function is determined; in the case of the first order, this distribution function is calculated as: In order for this scheme to be stable, it is necessary to fulfill the conditions on the basis of which the time step in the task is determined: To solve the relaxation equation, an explicit first-order difference scheme is used: The collision integral is calculated using the projection method described in more detail in [20]; the numerical integral obtained by this method is conservative in momentum, energy, and matter, and also vanishes on the Maxwell distribution. To increase the accuracy of calculating the collision integral, the nodes of the velocity grid for calculation are selected using the Korobov grid method [21].
The difference scheme is supplemented by initial and boundary conditions, the initial condition is the Maxwell distribution with a given pressure P 0 and temperature T 0 : Several types of boundary conditions are used in the work: diffuse reflection from the wall with a given temperature, constant gas flow, constant pressure condition, and specular reflection. To set the boundary conditions, it is necessary to determine the distribution function for molecules falling into the counting region through the boundary, i.e. for molecules with momenta (p, n) > 0, where n is the normal to the boundary directed inward into the counting region: where h(t) is the coefficient determined according to the type of boundary condition. In the case of the boundary condition of diffuse reflection, the coefficient h(t) is determined from the equality of gas flows to and from the wall, similarly, for the boundary condition of constant flow, the coefficient is determined based on the difference between the flow to and from the wall, determined by the specified value of the boundary flow, for the boundary condition of constant pressure coefficient is determined from normalization to a given pressure. In the case of a mirror boundary condition, the additionally determined distribution function has the form:

Programm
To implement the method of numerical solution of the Boltzmann equation, a software package [22] in C++ was developed that allows calculations for various initial and boundary conditions. Using the GMSH software package, the task geometry is set, the space is divided into cells of the required shape, and the cells are distributed into a given number of parallel nodes. To exchange data between parallel nodes, the openmpi library is used.

Algorithms
After setting the geometry of the problem, initial and boundary conditions, the spatial grid is divided between the specified number of parallel nodes. The separation process occurs using the GMSH package with the built-in Recursive or K-way algorithm, depending on the number of parallel nodes. For the counting algorithm to work correctly, fake cells are added to the spatial grid of each node, containing values of the distribution function of the cell of the neighboring node. After each data calculation operation, the distribution function values in the fake cells are synchronized. Since the fake cells of one process are neighboring the cells of the another process for which original cells has fake counterparts, synchronization is performed in two steps: • sending the values of the distribution function of all fake cells of the node to the neighboring node connected to the fake cell; • getting distribution function values for all fake cells.
A node with number 0 is used to input and output task parameters; this node is responsible for loading a spatial grid, creating dummy cells, and distributing it among other nodes. To obtain the final calculation results, macro parameters (pressure, temperature, etc.) in each cell of each node are collected on the main node and stored in a file.
6. Calculations 6.1. Calculations in the space of the IEG By the method described above, the problem of modeling the behavior of radionuclides in the space of the IEG was solved. Since the IEG has axial symmetry, to simplify the calculation, it was considered in a two-dimensional form, shown in Figure 1. The boundary conditions on the collector and emitter specify a diffuse reflection with a temperature of 900 K and 1500 K, respectively. The boundary condition at the left and right ends of the IEG is diffuse reflection  Figure 1. Two-dimensional view of a multi-element EGC: 1 -collector package; 2 -emitters; 3 -the flow of gases through the gas exhaust system with a gradient temperature and the constant pressure condition, respectively, for Cs the pressure is set to 150 Pa, for Kr and Xe the boundary condition for vacuum is set -pressure 0Pa. The initial condition is the Maxwell distribution with a pressure of 150 Pa for Cs and a gradient temperature, Kr and Xe are absent in the counting region at the initial moment. Since the estimated concentrations of Kr and Xe are small compared with the concentration of Cs, when calculating the collision integral, only collisions were considered: Cs -Cs, Cs -Kr and Cs -Xe, collisions of radionuclides with each other can be neglected. The calculation is made until the state of thermodynamic equilibrium is established. Serial calculations were carried out for different Kr and Xe flows in the IEG, and data were obtained on the distributions of pressures, temperatures, cesium and radionuclide flows in the gap cavity [23]. Figure 2 shows a plot of the flow of Kr and Xe at the right end of the IEG from the initial flow from the exhaust duct. At the exit from the IEG, the expected final flow qk was set. Using the approximation method, i.e. to minimize the squared error, the dependence corresponding to the obtained data is determined, from this dependence the value of the initial gas flow in the IEG is determined for a given expected output stream. The approximation was performed using the numpy package. Estimated in and expected out flows are shown in table 1. To estimate the partial pressures of Ba and Sr, calculations were performed with the initial flux determined by the approximation method; in the calculations, the beta decay of Kr and Xe to Sr and Ba was considered. In calculating the collision integral, only collisions were considered: Cs -Cs, Cs -Kr, and Cs -Xe. A total of 7 gases were considered: Cs, Kr, Xe, Rb, Sr, Cs (138), Ba. Figure 3 shows the pressure distribution graphs of Kr, Xe and Sr, Ba with respect to the length of the IEG. The gas dynamics calculation in the FSG was carried out using a three-dimensional grid consisting of prismatic cells, the schematic structure of the counting area is shown in Figure 4. Since the FSG has axial symmetry, to reduce the number of grid cells, a quarter of the cylinder is considered, mirror conditions are established on the symmetrical sides by mirror reflection. The initial condition was the Maxwell distribution at a pressure of 150 Pa and a temperature in the range 320-600 K for Cs, while Kr and Xe were initially absent in the counting region. A constant pressure Cs of 150 Pa and a flow of Xe and Kr were established as the boundary condition at the inlet of the FCG. At the exit from the FSG, the boundary condition of zero pressure -vacuum is established for all gases. In the region of evaporation and condensation, for Cs, the condition of a constant pressure of saturated vapors of 150 Pa is established at a temperature of 600 K, for Xe and Kr it is assumed that during the counting time they do not have time to dissolve in liquid cesium, therefore, all the stream entering the liquid cesium is transferred from the condensation region to evaporation area. On the walls of the FCG, the boundary condition sets the diffuse reflection with a temperature ranging from 600 K at the entrance to 320 K at the exit of the FCG. Due to the difference in estimated pressures of cesium and other gases, only collisions of pairs of molecules are considered: Cs -Cs, Cs -Kr and Cs -Xe. Serial calculations were performed for different Kr and Xe fluxes in the FCG, and data were obtained on the distributions of pressures, temperatures, cesium and radionuclide fluxes in the channel space. Figure 5 shows a graph of the dependence of the flow of Kr and Xe in the region of gas extraction -q k on the initial flow from the IEG -q 0 . Using the approximation method, gas flows were obtained at the entrance to the FCG for experimental output flows, determined from the data of the -spectrometric system. The received input and output streams are shown in table 2.

Conclusion
In the course of the work, the application of the method of numerical solution of the Boltzmann equation to the calculation of the behavior of the gaseous mixture in the cavity of the IEG and FSG was considered, a program was developed, and the distribution of radionuclide fluxes and