Nodal SN-method for HEX-Z geometry

JID: NUCET

ARTICLE IN PRESS

[m5+;May 6, 2016;9:53]

Available online at www.sciencedirect.com

Nuclear Energy and Technology 000 (2016) 1–4 www.elsevier.com/locate/nucet

Nodal SN-method for HEX-Z geometry ✩ V.P. Bereznev∗ Nuclear Safety Institute of the Russian Academy of Sciences, 52, B. Tulskaya Street, Moscow 115191, Russia Available online xxx

Abstract The problem of spatial approximation becomes very important in the solution of neutronics problems with coarse spatial grids, in particular, in the calculations of fuel assemblies of fast reactors (for instance, BN-800 and BN-1200 reactors) with computational cell in the form of hexagonal prism. “Weighted diamond difference” (WDD) schemes are the most widely used among the finite difference schemes for the neutron and gamma-ray transport equation solution. They are efficient from the viewpoint of ease of their implementation and associated CPU time expenditures. However, some drawbacks of these schemes are manifested when they are applied to solve the above described problems. Diamond difference scheme (DD) having second-order approximation (the best for this class of schemes) does not possess the properties of positivity and monotonicity. This is the reason why negative values and non-physical oscillations are often present in the solutions obtained. “Step” scheme (St), which is free from the disadvantages of the diamond difference scheme, has accuracy of only the first order. In connection with the need in high-accuracy calculations its use appears to be inefficient. There exist algorithms for correction of negative values, as well as adaptive (AWDD) schemes aimed both at the reduction of the level of oscillations and at the obtaining positive solutions. However, these algorithms negatively affect the order of approximation, and schemes of the first – second order of accuracy are discussed in such cases. Besides that, for adaptive schemes there exists the problem of correct selection of parameters of the scheme. The evident way to escape such situation with simultaneously enhancing quality and accuracy of the calculation is to select a fine mesh. In case of calculation of fuel assemblies of fast reactors spatial grid represents an arrangement of rectangular prisms with regular hexagons forming their bases (in such cases reference is made to HEX-Z-geometry). Therefore, hexagonal cells can be divided into rhomb-shaped cells (three rhombs per one hexagon; 12 rhombs per one hexagon, etc.). Diamond scheme is applied for the grids consisting of rhombs thus obtained. Because of the smaller cell size as compared with original cell size, the drawbacks inherent to this scheme will not be pronounced. Triangular grid can also be used. A different approach for the solution of the above indicated problem is to develop computational methods with enhanced order of accuracy without increasing the number of computational points. Nodal method is one of such methods. Expansion of unknown function inside the node (elementary volume with constant properties) in basis functions with subsequent calculation of expansion moments constitutes the basis of any nodal method. Nodal SN -method in HEX-Z geometry will be discussed in the present paper. Copyright © 2016, National Research Nuclear University MEPhI (Moscow Engineering Physics Institute). Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Keywords: Fast neutron reactor; Neutronics calculations; Nodal method; SN -approximation; Hexagonal geometry.

Nodal SN -method

✩ Russian text published: Izvestia Visshikh Uchebnikh Zavedeniy. Yadernaya Energetika (ISSN 0204-3327), 2015, n.3, pp. 56-62. ∗ Corresponding author. E-mail address: [email protected]. Peer-review under responsibility of National Research Nuclear University MEPhI (Moscow Engineering Physics Institute).

Let us examine stationary neutron transport process described by linear Boltzmann equation. After angular and energy discretization we obtain: g  g m ∇ϕm (x, y, z ) + (x, y, z )ϕmg (x, y, z ) t  G  M  χg  g g →g s = + ν f wm ϕmg (x, y, z ), (1) Ke f f m=1 g =1

http://dx.doi.org/10.1016/j.nucet.2016.03.004 2452-3038/Copyright © 2016, National Research Nuclear University MEPhI (Moscow Engineering Physics Institute). Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Please cite this article as: V.P. Bereznev, Nodal S N -method for HEX-Z geometry, Nuclear Energy and Technology (2016), http://dx.doi.org/10.1016/j.nucet.2016.03.004

ARTICLE IN PRESS

JID: NUCET

2

[m5+;May 6, 2016;9:53]

V.P. Bereznev / Nuclear Energy and Technology 000 (2016) 1–4

1 + μx



x

− x 2

[yb (t )Q(t ) − L(t )]

 t × exp − (x − t ) dt. μx

(3)

Now let us use polynomial expansion of neutron flux and the source as follows: (x) =

I 

I 

xin hi (x ), Q(x ) =

i=0

Qix hi (x),

i=0

 x 2 5 2 h(x) = 1, x, x − , . . . : yb (x)hi (x)h j (x) 12 − x 2

Fig. 1. Computational cell in HEX-Z-geometry.

where ϕmg (x, y, z ) is the neutron flux  density in point (x, y, z ) g in the direction m in group g; t (x, y, z ) is the total  macroscopic cross-section of the interaction; sg →g is the macroscopic cross-section of scattering from group g’ into group g; χ g is the fission neutron spectrum; Ke f f is the  effective neutron multiplication factor; ν gf is the number of fission neutrons emitted in one fission act; wm is the weight of the angular quadrature. Indices g and m will be hereinafter omitted by us when possible. Let us examine the computational cell with width across flats x and height zk as shown in Fig. 1 simulating one of the altitudinal sections of fuel assembly of fast reactor core:   x x × [−yb (x), yb (x)] × [zk−1/2 , zk+1/2 ] , D= − , 2 2

dx = δi j Dxj .

x − |x | yb (x ) = . √ 3

μ = μ cos α + η sin α,

Let us choose for the sake of certainty angular direction m > 0. Integration of Eq. (1) on variables y and z within the limits of the cell results in the following one-dimensional equation:

Thus, it is sufficient to make substitution (μ,η) → (μ ,η ) in expressions (3) and (4) for α u =π /3 and α v =2π /3. For variable z we have:

d yb (x)(x) + t yb (x)(x) = yb (x)Q(x) − L(x) dx For one-dimensional flux (x):

zk+1/2 yb (x) 1 (x) = dz ϕmg (x, y, z )dy. 2yb (x)zk zk−1/2 −yb (x) μ

(2)

Substitution of (4) in (3) allows obtaining expression for expansion moments and substitution of x = x/2 in (3) gives expression for flux xout on the outcoming plane. Applying similar technique we obtain expressions for functions (u) and (v) where variables u and v correspond to the coordinate axes shown in Fig. 1. When coordinate system is rotated by the angle α transfer equation in new coordinates (x´,y´,z´) is reduced to the following form:

  d  d  d μ + η + ξ + ϕ(x  , y , z ) = Q(x  , y , z ), dx dy dz t i.e. to Eq. (1) with following new direction cosines: η = −μ sin α + η cos α.

ξ d (z) + t (z) = Q(z) − L(z), z dz 2 (z) = √ 3x 2

L (x ) in Eq. (2) is the neutron leakage having radial component Lr (x) and axial component Lz (x) combined in (1) as follows [1]:

L(z) =

L (x ) = Lr (x) + yb (x )Lz (x )

(z) =

 μu u+ (x)−μv v+ (x) Lr (x) =

√ , 3 μv v+ (x)√ −μu u− (x) , 3

(4)

2 3x I 



x 2

− x 2



x <0

;

Lz (x) = ξ [z+ (x) − z− (x)]/zk . Taking into account the boundary condition (x/2 ) = xin we obtain the followig solution of Eq. (2):    t x x yb (x)(x) = √ exp − + x xin μx 2 2 3

yb (x)

d x ∫ ϕmg (x, y, z )d y, −yb (x)

μα [α+ (z) − α− (z)],

α={x,u,v}

zi fi (z ), Q(z ) =

i=0

x>0

(5)

I 

Qiz fi (z),

(6)

i=0

 zk+1/2 1 f (z) = 1, z, z2 − , . . . : fi (z) f j (z) dz = δi j Dzj . 12 zk−1/2 As the result we obtain the set of linear equations for neutron flux on the outcoming planes as follows: αout = a0α αin +

I  j=0

α = {x, u, v, z},

bα0 j L αj + c00 Qav +

I 

c0α j Qαj ,

j=1

(7)

Please cite this article as: V.P. Bereznev, Nodal S N -method for HEX-Z geometry, Nuclear Energy and Technology (2016), http://dx.doi.org/10.1016/j.nucet.2016.03.004

ARTICLE IN PRESS

JID: NUCET

[m5+;May 6, 2016;9:53]

V.P. Bereznev / Nuclear Energy and Technology 000 (2016) 1–4

3

Fig. 2. Map of the model of KNK-II reactor model [2].

where a0α , bα0 j , c0α j are the nodal coefficients and the following parity is taken into account:

Table 1 Results of calculation of Keff for two-dimensional problem.

Q0α = Qav ,

Calculation code (number

CR_in

of computation points)

Keff

K/K, %

Keff

K/K, %

where Qav is the average value of the neutron source inside the computational cell. Solution of the above set of equations is following:

MOCUM (90804) [3] SPANDOM (169) [4] Diamond scheme (169) Nodal scheme (169) TWOHEX-96 (16224)

1.00931 1.01055 0.99691 1.01009 1.00941

−0.01 0.11 −1.25 0.07 –

1.30866 1.30833 1.30620 1.30998 1.30945

−0.06 −0.09 −0.25 0.04 –

˜ xout ,  ˜ uout ,  ˜ vout ,  ˜ zout )T . out = (

CR_out

In this case the simplest algorithm for correction of negative values can be applied:

Table 2 Results of calculation of Keff for three-dimensional problem.

˜ αout ), αout = max (0, 

Calculation code (number

CR_in

of computation points)

Keff

K/K, %

Keff

K/K, %

HEXNOD [2] Diamond scheme Nodal scheme GMVP

1.0889 1.0911 1.0992 1.0955

−0.61 −0.4 +0.34 –

0.9783 0.9754 0.9876 −0.9839

−0.57 −0.87 +0.37 –

because balance equation is not used for calculation of fluxes on the outcoming planes and neutron balance will be taken into account during the phase of calculation of average flux av in the computational cell as follows: 2 3x +



|μα |(αout − αin ) +

α={x,u,v}



ξ (z − zin ) z k out

av = Qav .

Calculation of 2D- and 3D-models of KNK-II reactor (8)

t

Expansion moments of one-dimensional fluxes are calculated in accordance with the following expressions: αi = aiα αin +

I 

bαi j L αj + cαj0 Qav +

j=0

α = {x, u, v, z},

j = 1, . . . , I .

I 

ciαj Qαj ,

j=1

(9)

It is not necessary to calculate zero expansion moments because they are the same as the average value of flux in the computational cell: αi = av .

CR_out

Model of KNK-II reactor was chosen for testing (model No. 4 from [2]). This is the four-group conventionally critical problem. Two-dimensional model corresponding to the central altitudinal layer of the reactor core is shown in Fig. 2. Total number of fuel assemblies is equal to 169. Flat-to-flat dimension is equal to 12.99 cm. Two series of calculations were performed for the twodimensional problem (Table 1) corresponding to the presence (CR_in) and absence (CR_out) of control rods in the reactor core. Accuracy of convergence amounted for Keff to 10−6 . Calculation performed using TWOHEX-96 code [3] was taken as the reference. Results obtained using the nodal method are in good agreement with results obtained using other computer codes.

Please cite this article as: V.P. Bereznev, Nodal S N -method for HEX-Z geometry, Nuclear Energy and Technology (2016), http://dx.doi.org/10.1016/j.nucet.2016.03.004

JID: NUCET

4

ARTICLE IN PRESS

[m5+;May 6, 2016;9:53]

V.P. Bereznev / Nuclear Energy and Technology 000 (2016) 1–4

Fig. 3. Radial distribution of neutron field in the KNK-II model.

Calculations were also performed for three-dimensional model having 50 altitudinal layers. In this case state (CR_half) is added in which rods are partially submerged in the core. Calculation based on the use of Monte-Carlo method and the GMVP code [2] was chosen as the reference. As it is evident from Table 2 the nodal scheme outperforms the diamond scheme in terms of accuracy, especially in the case when heterogeneity is present. Radial distribution of neutron field presented in Fig. 3 demonstrates that the nodal scheme produces solutions without oscillations.

reactor performed in this study demonstrate enhanced order of accuracy and qualitatively correct distribution of neutron field when the nodal SN -method is used. References [1] H. Ikeda, T. Takeda, J. Nucl. Sci. Technol. 31 (1994) 497–509. [2] H. Ikeda, T. Takeda, 3-D Neutron Transport Benchmarks. Department of Nuclear Engineering, Osaka University, Japan, 1991 NEACRPL-330. [3] X. Yang, N. Satvat, Ann. Nucl. Energy 46 (2012) 20–28. [4] H. Tae, Z. Nam, Ann. Nucl. Energy 23 (1996) 133–143.

Conclusions Nodal SN -method was developed and applied for neutronics calculations. Calculations of 2D- and 3D- models of KNK-II

Please cite this article as: V.P. Bereznev, Nodal S N -method for HEX-Z geometry, Nuclear Energy and Technology (2016), http://dx.doi.org/10.1016/j.nucet.2016.03.004