Application of diffusion approximation in the calculations of reactor with cavities

The importance of calculation of radiation fields inside in-reactor cavities is associated with the necessity to simulate the emergency modes in fast breeder reactors (FBR), as well as reactor states with different coolant levels in special dedicated channels of passive feedback devices in lead-cooled fast reactors (LFR) of BREST type or in sodium cavities in sodium-cooled fast reactors (SFR). The Last Flight (LF) method (Bell and Glesston 1974, Davison 1960, DOORS3.2 1988, Mynatt et. al. 1969, Rhoades and Childs 1988, Rhoades and Sipmson 1997, SCALE 2009, Voloschenko et. al. 2012), or the method of the unscattered component is widely known and is commonly used in computer codes based on the method of spherical harmonics for obtaining solution in a gas medium at a certain distance from the calculated volume domain (DORT (Rhoades and Childs 1988), TORT (Rhoades and Sipmson 1997) and others (SCALE 2009)). The practice of its application (DOORS3.2 1988) demonstrated that acceptable results are achievable at considerable distances from the surface separating dense and gas media (more than two meters). Obtaining high-quality solution is not guaranteed for cavities within the calculation area. In addition, it is desirable to implement the cavities calculation methodology within the framework of the approximations used in reactor calculations introducing certain specific features. In particular isotropy of the neutron flux density and the necessity of forced introduction of a “conditional” calculation cell on the boundary surface of the void cavity are assumed in the diffusion approximation. If the LF method is oriented on the connection of the source point with the detection point, then it is necessary to determine in the calculation of neutron field in the cavities the neutrons escaping the surface area of the source and neutrons reaching a certain surface area of the cavity. In order to solve the problem, the authors suggested using the approximate solution presented in the paper. Thus, an algorithm for calculation of in-reactor cavities using the diffusion approximation was developed and implemented by the authors.


Introduction
Neutron transport in steady-state conditions in the fast reactor core can be described, for instance, in multi-group diffusion approximation as follows: where Ф g (r) is the neutron flux density; D g (r) is the diffusion factor; Σ g r (r) is the removal cross-section; Σ l→g (r) is the cross-section of neutron transfer from group l to group g; χ g is the fission neutron spectrum; νΣ g f (r) is the neutron multiplication cross-section; Q(r) is the intensity of the external neutron source.
The main advantage of the diffusion approximation for practical implementation in computer codes is the high speed with which the solution is obtained and low requirements on the availability of computational resources.Nevertheless, certain applicability limits exist posed by the presence of strongly absorbing medium and pronounced spatial heterogeneity.
In particular, diffusion approximation is not applicable in the presence of gaseous (vacuous) medium.Methods for calculating radiation in such media can be based on the solution of integral neutron transport equation.If the purpose is to define radiation field outside the calculation area, then the last flight method (LF) is widely applied.However, in practice this method is not easily applicable for calculating internal cavities, since uncertainty of the method is high at the distances close to the surface.In such case volume sources used in LF-method can be replaced with surface sources with subsequent numerical integration over the surface separating dense and gaseous media.
Thus, the purpose of the present study is to create combined algorithm with associated calculation module converging the diffusion approximation producing fast solution with semi-analytical solution of integral equation inside in-reactor cavities.

Algorithm for calculating unscattered component
Linear integrodifferential neutron transport equation can be written in the following integral form: (1) Equation ( 1) means that neutron flux φ(r, Ω, E, t) in point r is caused by neutrons which appeared in all points r-s′Ω with direction Ω and energy E for all positive s′.Exponential multiplier here is the attenuation factor characterizing the reduction of neutron flux when s = 0 is reached, and Σ(r-s′′Ω, E) is the total macroscopic cross-section of neutron interaction with medium.Integration over s′ can be achieved only within the bounda-ries of the area under examination if incoming neutron flux is absent.
We suggest spreading such representation of the solution only on "vacuous" calculation cells filled with gaseous medium where density of the matter is low to such extent that cross-sections of neutron interactions with medium are very small, so that diffusion factor exceeds by several orders of magnitude its values in a typical reactor medium, and not allowing using conventional diffusion approximation in the medium in question.
As it was indicated in (Davison 1960), when the non-zero flux of incoming neutrons is present, it can be replaced with equivalent surface source.Such transition is possible due to the calculation of neutron flux density on the boundary surface between dense and gaseous media.where integration is performed over the surface separating dense and gaseous media; exp(-Σ g |r -r′|)/(4π|r -r′| 2 ) is the probability for neutron from r to reach r′ without collisions; φ g (r′)dS′ is the field from the surface element dS′; φ g (r) is the field in r due to the field from surface S′.Multi-group energy approximation to which index g corresponds is used in expression (2) and hereafter.
Let us examine transmission of radiation using two-dimensional model of hexagonal calculation cell conventionally applied in the FBR calculations as represented in Figure 1 where integration is performed along the plane AB; φ g (r) is the field in point r taking into account radiation coming from only the AB plane.
Average value of field on the "receiving" plane ab taking into account (3) and rearrangement of the order of integration is as follows: In such case the value is the contribution in the field on the plane ab made only from the field in point r′ located on the plane AB.
Assuming that field on the plane AB is the same for all its points, i.e. that φ g (r′) = φ g AB , and y -y′=const=Δy in accordance with Fig. 1, we obtain: The obtained integral can be calculated approximately and, in case of gaseous medium (Σ g ~10 -4 cm -1 ), Taylor expansion of the integration term can be used with ( 6) accepting the following form: In the simplest case for x -x′ = const = Δx expression (6) attains the following form: exp , 2 where s is the distance between the points of the "emitting" and "accepting" surfaces depending in the general case on the angle at which the point on the accepting surface is visible from the point of emitting surface, i.e. s = s(α).
Since the problem of radiation transmission in the hexagonal cell (see Fig. 1) is symmetrical relative to vertical axis passing through the point S, then the solution (8) can be determined only for the length AS, after which the obtained value must be doubled for the whole length AB: exp .
Neutron flux density on the boundary separating the media can be determined from the integrodifferential neutron transport equation written for the boundary surface, i.e. for the zero-measure volume during integration of which all members associated with cell volume will disappear and only the flow members, integral for which will be converted from volume to surface ones, will remain.In such case we obtain equality of flows on the surface separating media exp 4 exp , 4 where φ g i is the neutron flux density on i-th surface of the cavity from which neutron can reach the analyzed surface s; α is , α si are the apertures of angles under which surface s is seen from the surface i and vice versa; Δ is is the distance between the surfaces i and s; s j is the distance in the dense cell where diffusion approximation works with diffusion factor D j between the center-of-mass point with neutron flux density φ g i and the boundary surface between the media s; S i , S s are the surface areas of calculated cells i and j with shared gaseous medium; φ g s is the neutron flux density on the surface s separating media in the calculated cell j.
Since integration strongly depends on the cell geometry, then let us describe them.

Types of calculation cells with voids
There are cells with external steel cladding and central cells with coolant and some structures (end fittings of fuel rods in reactors of BN type (Saraev et. al. 2010, Oshkanov et. al. 2010) or technological pipes, as, for instance, in Figure 2) homogenously distributed in it.
While the above mentioned homogenous medium fills the channel in the fuel assembly of BN-type reactor in normal operational conditions starting from the cells containing absorber to the top of reactor core and is substituted with gaseous medium in the case of boiling coolant, then in normal operational conditions of BREST-type reactors (Dragunov et. al. 2012) the above homogenous medium fills the positive feedback devices to the reactor core top and gives the place to the gaseous medium (argon) when pressure in the lower header of the reactor primary cooling loop drops with the forward stroke of the gas of 70 cm from the top of the core and lower.
Corresponding model of calculation cells and nodes is represented in Figure 2: for seven calculation cells there are seven main nodes located in the centers of mass of the cells and six auxiliary nodes arranged along the internal boundaries of the cells.If auxiliary nodes along the height are added one for each internal boundary of internal hexagonal calculation cells, then we obtain seven auxiliary calculation nodes for each altitudinal row of the channel.In case of complete coolant drainage from calculation cells of the central channel and their filling with inter gas, diffusion factor in such medium varies from unity to ~100 for positive feedback devices and to ~10 4 for reactor facilities of BN type.Conventional diffusion approximation does not provide for the obtaining of acceptable result.Auxiliary calculation nodes become determining, the problem of obtaining average diffusion factor on the boundaries of calculation cells is removed and solution of integral neutron transport equation in gaseous medium is ensured.
Possible model of the positive feedback device includes external steel cladding and central steel technological pipe between which coolant is located filling the channel to the reactor core top in normal operational conditions and giving place to gaseous medium (argon) when pressure drops in the lower header of the first reactor cooling loop with forward gas stroke of 70 cm from the reactor core top and lower (Fig. 3).
Model of calculation cells and nodes in the positive feedback device channel is also represented in Figure 3: for 13 calculation cells there are 13 main nodes arranged in the centers of mass of the cells, and 12 auxiliary nodes located on internal boundaries of the cells.If auxiliary nodes which are arranged along the height one for each internal altitudinal boundary of internal trapezoidal calculation cells are added to them, then we obtain 18auxiliary calculation nodes for each altitudinal row of the channel.
Such layout of arrangement of calculation nodes allows using conventional procedure in diffusion approximation with filling intermediate calculation cells with coolant when auxiliary calculation nodes are actually playing supporting role.However, when coolant is drained from auxiliary calculation cells calculation nodes on the surface separating media play determining role in the solution of neutron transport equations in gaseous medium.
Thus, concordance is achieved for the indicated layout between the solution of diffusion equation in dense media and solution of integral equation in gaseous medium.
Specific models of radiation in calculation cells filled with gaseous medium are examined below.

Algorithm of calculation of void in hexagonal cell
Let us examine possible situations in the propagation of radiation in hexagonal cell with internal hexagonal void.
When hexagonal cell (Fig. 4) is drained, neutrons are drifting from surface 1 to surfaces 2 -6, and arrival of neutrons from these surfaces will occur correspondingly.
Values of angles for points А, S, and В are represented as аСn, where C is the name of the calculation cell, n is the number of the receiving surface.Lengths of sections are denominated similarly to (xCn).
When equation ( 9) is used for calculations, the initial section АВ must be divided for improving accuracy, for instance, into 10 sections with dimensions equal to АВ/10 each, which significantly reduces variations of the aspect angle and average distance from the points on the section emitting radiation to the receiving surface.Besides the above we assume linear dependence of flight distance on the angle under which the point of reception of radiation represented in Figure 5 is visible from the point of emission of radiation (х 0 , х 1 , х 2 are the distances from the source to the receiving side, х 0 is the shortest distance; α 1 , α 2 are the aspect angles under which receiving sections ab and bc are visible from the point of the source).

Layout of propagation of radiation in the void between the external and internal dense hexagonal structures
Let us examine possible situations with propagation of radiation in the hexagonal cell with gas void located between the external and internal dense media represented in hexagonal form.
First let us examine propagation of radiation from the internal dense hexagonal structure.For surface 7 with flat-to-flat dimension equal to h c in hexagonal cell with internal hexagonal cell drift of neutrons will occur to the surfaces 1, 2, 6 both from the angle, as well as from all other, for instance, central points of surface 7, which is shown in Fig. 6.
Let us analyze backward radiation from the external hexagonal structure to the internal one through gaseous medium.For this purpose, we will examine propagation of radiation for points of one of the side planes, for instance, the first one, depending on the dimensions of the internal hexagonal cell with flat-to-flat size h c relative to the dimensions of the external hexagonal cell with flat-to-flat dimension equal to g. Varying the ratio of parameter h c to g we obtain variation of visibility of the sides of the cell under examination from points on the side plane 1 (AB).These variations are observed on the boundaries of sections with the indicated values of the ratio: 2/3g ≤ h c ; g/2 ≤ h c < 2/3g; g/3 ≤ h c < g/2; g/4 ≤ h c < g/3; g/6 ≤ h c < g/4; h c < g/6.
Propagation of radiation from points of the side plane AB is presented in Figure 7 for h c < g/6.In this case point D becomes visible from point P located somewhere between points G и V with AV > AP and half-line РD touching the angle of internal hexagonal element.Thus, side plane 4 becomes visible from points of the section PS in the area of interface with side plane 3, and side plane 8 becomes visible from points of section GS 8.
Point А in Figure 7 is the extreme left point of the side plane 1; S is its midpoint; V; G, P are the intermediate points between A and S. For the case under examination side planes 2, 3, 4a, 5, 6, 7 and 12 are visible from points of section AG; side planes 2, 3, 4а, 5, 6, 7, 12 and 8 are visible from points of section GP; side planes 2, 3, 4a and 4b, 5, 6, 7, 8 and 12 are visible from points of section PV; side planes 2, 3, 4a and 4b, 5, 6, 7, 8 and 12 are visible from points of section VS.

Preliminary results of application of the algorithm of calculation of reactor with voids
Results of comparison of design-basis drainage of 18 positive feedback devices in the version of design of BREST reactor facility obtained using Monte-Carlo method (Alekseev et. al. 2016, Gurevich et. al. 2016) and the above described algorithm are presented in Table 1.
Calculation time for model of BREST-type reactor facility with model of positive feedback device presented   in Fig. 3 did not exceed 25 seconds on one computer processor core with accuracy equal to 10 -6 for neutron flux density of multi-group diffusion approximation.

Conclusion
Description is given of the algorithm for calculation in diffusion approximation of FBR in the presence of calculation cells containing voids.The algorithm for solution of neutron transport equation in gaseous medium is based on the semi-analytical methodology of solution of Peierls integral equation in the assumption that calculation cells with gaseous medium are surrounded by calculation cells in which diffusion approximation is used, which allows keeping high speed of obtaining the solution with high enough accuracy of this solution.No information about earlier application of such algorithm in international calculation practice is known to authors.
Preliminary results of application of the suggested algorithm are of high enough quality with high speed of obtaining the solution.As it was originally expected, direct application of diffusion approximation produces solution with accuracy which is not acceptable.
The developed algorithm for calculation in diffusion approximation of fast breeder reactor with presence of calculation cells with voids can be used for estimation of sodium void reactivity effect (Poplavsky et. al. 2010) for any reactor facility equipped with FBR, as well as for estimation of behavior of neutron field in FBR with disintegrating core in the analysis of postulated accidents for assessment of safety of the specific reactor facility.

Figure 1 .
Figure 1.Model of transmission of radiation in hexagonal calculation cell with internal hexagonal structure (S is the center point of the length AB).

Figure 2 .
Figure 2. Model of calculation cell with drainable internal hexagonal cell and arrangement of calculation nodes.

Figure 3 .
Figure 3. Model of calculation cell of the positive feedback device of BREST reactor facility and arrangement of calculation nodes.

Figure 4 .
Figure 4. Propagation of radiation from points А, S and B of the hexagonal cell (AS = AB/2).

Figure 5 .
Figure 5. Dependence of flight distance on the aspect angle of the receiving section.

Figure 6 .
Figure 6.Propagation of radiation from the points of section ab on the side 7.

Figure 7 .
Figure 7. Propagation of radiation from points A, G, P, V and S.

Table 1 .
Deviation of efficiency of drainage of the positive feedback device to the height of 70 cm from the reactor core top from the designed value, %.