2. Reactor core calculation overview

1. Introduction

6 New Trends in Nuclear Science

concurrently.

prompt results.

The mathematical models representing the nuclear reactor physics are based mainly on two theoretical areas: neutron transport theory and neutron diffusion theory, where it is necessary to remark that neutron diffusion theory is really a simplification of the neutron transport theory.

Numerical methods are used to solve the partial differential equations representing the nuclear reactor physics, and these methods are derived from discretization techniques. For numerical solutions in any scientific area, computational tools have been developed including software and hardware. In the past, the former computer processing was the sequential execution of computer commands, meaning to say that program tasks are carried out one after one. Modern computational tools have been developed for parallel processing, executing several tasks

The computing branch dealing with the system architecture and appropriate software related to the simultaneous execution of computer instructions and applications is known as parallel computing science. Former developments in parallel computing were made in the late 1950s, following the construction of supercomputers throughout the 1960s and 1970s. Nowadays, clusters are

Since the late 1950s, the performance of safety analyses was essential in the nuclear industry, in research reactors, but mainly safety analyses of nuclear power plants for commercial purposes. Scientific computing calculations were vital to these safety analyses, but with important limitations in computer/computing capabilities. At the beginning, the objective was to give a solution to partial differential equation models based on neutron diffusion or neutron transport with technology and methods available in those years. Numerical techniques were used first with finite differences and finite element approaches, and gradually up to now, with nodal finite element methods (NFEMs). Despite the numerical method employed, the computer code user faces the problem of solving extremely large algebraic systems challenging hardware/ software capabilities. Generation of results for any reactor simulation in considerable short

Recent developments of high-performance computer equipment and software have made the use of supercomputing in many scientific areas possible. The appropriate selection of parallel computing software, like newly developed linear algebra libraries, to be used in a specific project may result in a suitable platform to simulate nuclear reactor states with relatively

Throughout the world, several research projects in the last decade have been developed with the main objective of making full tridimensional (3D) coupling simulations of nuclear reactor cores, leaving aside the obsolescence of the point kinetics theory. Most of the modern nuclear reactor simulators are based on neutron transport theory, or on neutron diffusion theory, to obtain detailed 3D results. As light water is used for cooling/moderating light water reactors (LWRs), a comprehensive analysis of the reactor core physics must include thermal-hydraulic phenomena, so that modern simulations perform reactor calculations with thermal-hydraulic

the workhorse of scientific computing and are the dominant architecture in data centers.

times is a desirable achievement for computer code users [1].

feedback coupled with neutron kinetics calculations.

Although there has been growing interest in the transport-based core neutronics analysis methods for a more accurate calculation with high-performance computers, it is yet impractical to apply them in the real core design activities because their performance is not so practical on ordinary desktop or server computing machines. For this reason, most of the neutronics codes for reactor core calculations are still subject to the two-step calculation procedure, which consists of (1) homogenized group neutron parameters generation and (2) neutron diffusion core calculation.

In the core calculation steps that are the main concern of this work, nodal codes based on the diffusion theory have been used to determine the neutron multiplication factor and the corresponding core neutron flux (or power) distribution. Practically, almost all nuclear reactor simulation codes employ the two-group approach involving only fast and thermal neutron energy groups for the applications to light water reactors (LWR). However, numerical calculations with the two-group structure are not appropriate in the analysis of cores loaded with mixed oxide fuels or analysis of fast breeder reactors, since the neutron spectrum is influenced more by the core environment, requiring much more energy groups than only two groups.

As settled in Ref. [2], even using a high-performance computer, a direct core calculation with several tens of thousands of fuel pins is difficult to perform in its heterogeneous geometry model form, using fine groups of a prepared reactor cross-section library. The Monte Carlo method can handle such a core calculation (see also the Serpent code), but it is not easy to obtain enough accuracy for a local calculation or small reactivity because of accompanying statistical errors, besides the large calculation times. Instead of using neutron transport computer codes, the nuclear design calculation is performed in two steps: (1) lattice calculation in a two-dimensional infinite arrangement of fuel rods or assemblies for the generation of homogenized lattices jointly with their corresponding homogenized cross-sections and (2) core calculation in a three-dimensional whole core, with a neutron diffusion code using the information of the previous step.

As shown in Figure 1 [2], the lattice calculation prepares few-group homogenized cross sections which maintain the energy dependence (neutron spectrum) of nuclear reactions, and these reduce the core calculation cost in terms of time and memory. The final core design

<sup>1</sup> This work was performed under the auspices of the financial support from the National Strategic Project No. 212602 (AZTLAN Platform) as part of the Sectorial Fund for Energetic Sustainability CONACYT—SENER, Mexico.

Figure 1. Typical lattice calculation process flow for light water reactors [2].

parameters are not concerned with continuous energy dependence, but spatial dependence, such as power distribution, is important to avoid high local neutron fluxes or high absorbing materials causing significant neutron flux gradients, mainly when safety analyses are performed upon the final proposed core designs.

In the core calculations with space-dependent data (cross sections and neutron flux), the effective cross sections are processed, with a little degradation in the accuracy as possible, by using the results from the multi-group lattice calculation. Lattice code calculation and codes are not discussed here.

3. Neutron diffusion theory and nodal methods

Figure 2. Homogenization and group collapsing of cross sections [2].

1 vg ∂ ∂t ϕ<sup>g</sup> r !; t

∂ ∂t Ci r !; t � � <sup>¼</sup> <sup>β</sup><sup>i</sup>

� � <sup>¼</sup> <sup>∇</sup>•D<sup>g</sup>∇ϕ<sup>g</sup> <sup>r</sup>

X G

g¼1

are described in [5].

<sup>þ</sup> <sup>1</sup> � <sup>β</sup> � �χ<sup>g</sup>

ν<sup>g</sup> r !; t � �Σ<sup>g</sup>

3.1. Multi-group time-dependent neutron diffusion equations

!; t � � � <sup>Σ</sup><sup>g</sup>

X G

> <sup>f</sup> r !; t � �ϕ<sup>g</sup> <sup>r</sup>

g0 ¼1 νg0 r !; t � �Σ<sup>g</sup><sup>0</sup>

For G neutron energy groups and Ip delayed neutron precursor concentrations, the neutron diffusion kinetics equations are given by Eqs. (1) and (2) [5]. Although there has been a growing interest in the transport-based core neutronics analysis methods for more accurate calculation with high-performance computers, it is yet impractical to apply them in the real core design activities because their performance is not so practical on ordinary desktop or server computing machines. For this reason, most of the neutronics codes for reactor core calculations are still subject to the two-step calculation procedure, which consists of homoge-

> !; t � � <sup>þ</sup> <sup>X</sup>

> > r !; t � � þ<sup>X</sup>

!; t

G

Σg0 !<sup>g</sup> <sup>s</sup> r !; t � �ϕ<sup>g</sup><sup>0</sup>

Ip

i¼1 χg <sup>i</sup> λiCi r !; t

� �; i <sup>¼</sup> <sup>1</sup>,…, Ip; <sup>∀</sup> <sup>r</sup>

r !; t � �

!; t

� �; <sup>g</sup> <sup>¼</sup> <sup>1</sup>, …, G;

� �<sup>∈</sup> <sup>Ω</sup> � ð � <sup>0</sup>; <sup>T</sup> (2)

Nuclear Reactor Simulation

9

http://dx.doi.org/10.5772/intechopen.79723

(1)

g<sup>0</sup> ¼ 1 g<sup>0</sup> 6¼ g

nized group neutron parameter generation and neutron diffusion core calculation

<sup>R</sup> r !; t � �ϕ<sup>g</sup> <sup>r</sup>

> <sup>f</sup> r !; t � �ϕ<sup>g</sup><sup>0</sup>

In addition to boundary conditions for neutron fluxes, initial conditions must be satisfied by neutron fluxes and neutron precursor functions. Parameters involved in the above equations

!; t � � � <sup>λ</sup>iCi <sup>r</sup>

There are two processes followed for lattice calculation. One is the homogenization to lessen the space-dependent information and the other is group-collapsing to reduce the energydependent information as shown in Figure 2. The fundamental idea of both methods is to preserve neutron reaction rate. The next step is to consider the conservation of reaction rate in the energy group G in the same manner as that in the homogenization.

The number of few groups depends on reactor type and computation code. Two or three groups are adopted for the NK- and TH-coupled core calculation of LWRs and much more groups (18, 33, etc.) are used for the core calculation of LMFRs (Liquid Metal Fast Reactors). Currently, revised methods exist for the improvement of cross-sections generation using computer codes dedicated to lattice calculation for few-groups approach, like in Ref. [3], where three topics are involved: (1) improved treatment of neutron-multiplying scattering reactions; (2) group constant generation in reflectors and other non-fissile regions, leading to the use of discontinuity factors in neutron diffusion codes; and (3) homogenization in leakage-corrected criticality spectrum, in which several leakage corrections are used to attain criticality, accounting for the non-physical infinite-lattice approximation. Another improvement was done in Monte Carlo codes [4], implementing reliable multi-group cross-sections calculations for collapsed flux spectrum. Ref. [4] focuses on calculating scattering cross sections, including the group-to-group scattering.

The following sections contain, as a matter of example, summarized explanations of the AZKIND nuclear reactor simulator in which the reactor physics is based on neutron diffusion theory.

Figure 2. Homogenization and group collapsing of cross sections [2].

parameters are not concerned with continuous energy dependence, but spatial dependence, such as power distribution, is important to avoid high local neutron fluxes or high absorbing materials causing significant neutron flux gradients, mainly when safety analyses are performed

In the core calculations with space-dependent data (cross sections and neutron flux), the effective cross sections are processed, with a little degradation in the accuracy as possible, by using the results from the multi-group lattice calculation. Lattice code calculation and codes

There are two processes followed for lattice calculation. One is the homogenization to lessen the space-dependent information and the other is group-collapsing to reduce the energydependent information as shown in Figure 2. The fundamental idea of both methods is to preserve neutron reaction rate. The next step is to consider the conservation of reaction rate in

The number of few groups depends on reactor type and computation code. Two or three groups are adopted for the NK- and TH-coupled core calculation of LWRs and much more groups (18, 33, etc.) are used for the core calculation of LMFRs (Liquid Metal Fast Reactors). Currently, revised methods exist for the improvement of cross-sections generation using computer codes dedicated to lattice calculation for few-groups approach, like in Ref. [3], where three topics are involved: (1) improved treatment of neutron-multiplying scattering reactions; (2) group constant generation in reflectors and other non-fissile regions, leading to the use of discontinuity factors in neutron diffusion codes; and (3) homogenization in leakage-corrected criticality spectrum, in which several leakage corrections are used to attain criticality, accounting for the non-physical infinite-lattice approximation. Another improvement was done in Monte Carlo codes [4], implementing reliable multi-group cross-sections calculations for collapsed flux spectrum. Ref. [4] focuses on calculating scattering cross sections, including the

The following sections contain, as a matter of example, summarized explanations of the AZKIND nuclear reactor simulator in which the reactor physics is based on neutron diffusion theory.

the energy group G in the same manner as that in the homogenization.

upon the final proposed core designs.

Figure 1. Typical lattice calculation process flow for light water reactors [2].

are not discussed here.

8 New Trends in Nuclear Science

group-to-group scattering.

### 3. Neutron diffusion theory and nodal methods

#### 3.1. Multi-group time-dependent neutron diffusion equations

For G neutron energy groups and Ip delayed neutron precursor concentrations, the neutron diffusion kinetics equations are given by Eqs. (1) and (2) [5]. Although there has been a growing interest in the transport-based core neutronics analysis methods for more accurate calculation with high-performance computers, it is yet impractical to apply them in the real core design activities because their performance is not so practical on ordinary desktop or server computing machines. For this reason, most of the neutronics codes for reactor core calculations are still subject to the two-step calculation procedure, which consists of homogenized group neutron parameter generation and neutron diffusion core calculation

1 vg ∂ ∂t ϕ<sup>g</sup> r !; t � � <sup>¼</sup> <sup>∇</sup>•D<sup>g</sup>∇ϕ<sup>g</sup> <sup>r</sup> !; t � � � <sup>Σ</sup><sup>g</sup> <sup>R</sup> r !; t � �ϕ<sup>g</sup> <sup>r</sup> !; t � � <sup>þ</sup> <sup>X</sup> G g<sup>0</sup> ¼ 1 g<sup>0</sup> 6¼ g Σg0 !<sup>g</sup> <sup>s</sup> r !; t � �ϕ<sup>g</sup><sup>0</sup> r !; t � � <sup>þ</sup> <sup>1</sup> � <sup>β</sup> � �χ<sup>g</sup> X G g0 ¼1 νg0 r !; t � �Σ<sup>g</sup><sup>0</sup> <sup>f</sup> r !; t � �ϕ<sup>g</sup><sup>0</sup> r !; t � � þ<sup>X</sup> Ip i¼1 χg <sup>i</sup> λiCi r !; t � �; <sup>g</sup> <sup>¼</sup> <sup>1</sup>, …, G; (1) ∂ ∂t Ci r !; t � � <sup>¼</sup> <sup>β</sup><sup>i</sup> X G g¼1 ν<sup>g</sup> r !; t � �Σ<sup>g</sup> <sup>f</sup> r !; t � �ϕ<sup>g</sup> <sup>r</sup> !; t � � � <sup>λ</sup>iCi <sup>r</sup> !; t � �; i <sup>¼</sup> <sup>1</sup>,…, Ip; <sup>∀</sup> <sup>r</sup> !; t � �<sup>∈</sup> <sup>Ω</sup> � ð � <sup>0</sup>; <sup>T</sup> (2)

In addition to boundary conditions for neutron fluxes, initial conditions must be satisfied by neutron fluxes and neutron precursor functions. Parameters involved in the above equations are described in [5].

#### 3.2. Spatial discretization

The spatial discretization of Eqs. (1) and (2) is strongly connected with the discretization of a nuclear reactor core of volume Ω. Representing the neutron flux and the precursor concentrations in terms of base functions defined over Ω, it is possible to write

$$\phi^{\mathfrak{g}}\left(\vec{r},t\right) \equiv \sum\_{k=1}^{N\_f} u\_k\left(\vec{r}\right) \phi\_k^{\mathfrak{g}}(t);\ \mathbf{g} = 1,...,G;\ \forall \left(\vec{r},t\right) \in \Omega \times (0,T];\tag{3}$$

3.3. NFE method in spatial discretization

[xi

,xi+1] � [yj

u<sup>00</sup>

u<sup>00</sup>

u<sup>00</sup>

where Plpqð Þ¼ x; y; z Plð Þx Ppð Þy Pqð Þz .

<sup>L</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> ¼ � <sup>1</sup>

<sup>N</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> ¼ � <sup>1</sup>

<sup>B</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> ¼ � <sup>1</sup>

u<sup>000</sup>

be written for the NFE method RTN-0 (Raviart-Thomas-Nédélec).

2

2

2

As fully explained in [6] and summarized in [1], a simple NFE element is characterized by the fact that for each node, the function unknowns to be determined are the (00) Legendre moment (average) of the unknown function over each face of the node and the (000) Legendre moment over the node volume. Figure 3(a) shows a physical domain Ω graphically represented after generating an xyz mesh. Figure 3(b) shows a cuboid-type node with directions through the faces: (x) Right, Left; (y) Near, Far; (z) Top, Bottom; and C for the average of the function over the node volume. Taking into consideration the general form to build up nodal schemes [7], the moments of a function (at edges and body) over a node like the one shown in Figure 3(b) can

In the NFE method RTN-0, the normalized zero-order Legendre polynomials defined over the unit cell Ωijk = [�1,+1] � [�1,+1] � [�1,+1] and correlated to each physical cell Ω<sup>e</sup> = Ωijk =

The matrix elements are quantified introducing the following nodal basis functions [7]:

ð Þ P100 � P200 ; <sup>u</sup><sup>00</sup>

ð Þ P010 � P020 ; <sup>u</sup><sup>00</sup>

ð Þ P001 � P002 ; <sup>u</sup><sup>00</sup>

Figure 3. Discretization of reactor volume Ω and a local node Ωe. (a) Domain Ω. (b) Physical local node Ωe.

<sup>C</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> <sup>¼</sup> P000 � P200 � P020 � P002;

,yj+1] � [zk,zk+1] are used to calculate the elements of the matrices in Eqs. (5) and (6).

<sup>R</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> ¼ þ

<sup>F</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> ¼ þ

<sup>T</sup> <sup>x</sup>; <sup>y</sup>; <sup>z</sup> ¼ þ

1 2

1 2

1 2 ð Þ P100 þ P200 ;

ð Þ P001 þ P002 ;

ð Þ P010 þ P020 ; (7)

Nuclear Reactor Simulation

11

http://dx.doi.org/10.5772/intechopen.79723

$$\mathbf{C}\_{i}\left(\overrightarrow{r},t\right) \equiv \sum\_{m=1}^{N\_{p}} \mathbf{v}\_{m}\left(\overrightarrow{r}\right) \mathbf{C}\_{i}^{m}(t); \ i = 1, \ldots, I\_{p}; \ \forall \left(\overrightarrow{r},t\right) \in \Omega \times (0,T]; \tag{4}$$

where Nf and Np are the number of unknowns to be determined for neutron flux and delayed neutron precursors, respectively. Substituting expressions (3) and (4) into (1) and (2), and applying the Galerkin process for spatial discretization, as described in [6], the resulting algebraic system of equations can be expressed in a matrix notation as follows:

$$\begin{aligned} \frac{1}{\mathbf{c}^{\mathfrak{g}}} \mathbf{M}\_{\uparrow} \frac{d}{dt} \phi^{\mathfrak{g}}(t) &= -\mathbf{K}^{\mathfrak{g}} \phi^{\mathfrak{g}}(t) - \sum\_{\mathbf{g}'=1}^{\mathbf{G}} \mathbf{S}^{\mathbf{g}' \rightarrow \mathbf{g}} \phi^{\mathfrak{g}'}(t) \\ &+ (1 - \mathfrak{f}) \chi^{\mathfrak{g}} \sum\_{\mathbf{g}'=1}^{\mathbf{G}} \mathbf{F}^{\mathbf{g} \mathbf{g}'}(t) \phi^{\mathfrak{g}'}(t) + \sum\_{i=1}^{l\_{\mathfrak{f}}} \mathbf{H}^{\mathfrak{g} i}(t) \mathbf{C}\_{i}(t), \ \mathbf{g} &= 1, \dots, \mathbf{G}; \\\\ \mathbf{d}\_{\mathbf{g}', \mathbf{g}'} &= \sum\_{\mathbf{g}', \mathbf{g}'}^{\mathbf{G}} \mathbf{H}^{\mathbf{g}', \mathbf{g}}(t) \phi^{\mathfrak{g}', \mathbf{g}'}(t), \ \mathbf{C}\_{i}(t), \ \mathbf{c}\_{i}(\pm 1, \pm \omega, \dots, \omega, \pm 1, \dots, \omega, \dots, \omega, \dots, \omega, \dots, \omega, \dots, \omega, \dots, \omega, \dots, \omega, \dots, \omega. \end{aligned} (5)$$

$$\mathbf{M}\_p \frac{d}{dt} \mathbf{C}\_i(t) = \sum\_{\mathbf{g'}=1}^{\bullet} \mathbf{P}^{\mathbf{g}'}(t) \boldsymbol{\phi}^{\mathbf{g'}}(t) - \lambda\_i \mathbf{M}\_p \mathbf{C}\_i(t), \ \mathrm{i} = 1, \dots, I\_p; \ \forall \left(\overrightarrow{r}, t\right) \in \Omega \times (0, T]; \tag{6}$$

where <sup>ϕ</sup><sup>g</sup>ðÞ¼ <sup>t</sup> <sup>ϕ</sup><sup>g</sup> <sup>1</sup> ð Þ<sup>t</sup> ;…; <sup>ϕ</sup><sup>g</sup> Nf ð Þt h i<sup>T</sup> and CiðÞ¼ <sup>t</sup> <sup>C</sup><sup>1</sup> <sup>i</sup> ð Þ<sup>t</sup> ; …;CNp <sup>i</sup> ð Þt h i<sup>T</sup> : Table 1 contains the expressions representing the calculation of each matrix coefficient.


Table 1. Matrix elements from the spatial discretization.

#### 3.3. NFE method in spatial discretization

3.2. Spatial discretization

10 New Trends in Nuclear Science

ϕ<sup>g</sup> r !; t � �

Ci r !; t � �

1 vg M<sup>f</sup> d

M<sup>p</sup> d

where <sup>ϕ</sup><sup>g</sup>ðÞ¼ <sup>t</sup> <sup>ϕ</sup><sup>g</sup>

Sg0

dt CiðÞ¼ <sup>t</sup> <sup>X</sup>

G

g0 ¼1 P<sup>i</sup>g<sup>0</sup> ð Þ<sup>t</sup> <sup>ϕ</sup><sup>g</sup><sup>0</sup>

<sup>1</sup> ð Þ<sup>t</sup> ;…; <sup>ϕ</sup><sup>g</sup>

Table 1. Matrix elements from the spatial discretization.

h i<sup>T</sup>

Nf ð Þt

sions representing the calculation of each matrix coefficient.

The spatial discretization of Eqs. (1) and (2) is strongly connected with the discretization of a nuclear reactor core of volume Ω. Representing the neutron flux and the precursor concentra-

<sup>k</sup> ð Þt ; g ¼ 1, …, G; ∀ r

<sup>i</sup> ð Þt ; i ¼ 1, …, Ip; ∀ r

ð Þþ <sup>t</sup> <sup>X</sup> Ip

<sup>i</sup> ð Þ<sup>t</sup> ; …;CNp

h i<sup>T</sup>

<sup>i</sup> ð Þt

ð Þ� t λiMpCið Þt , i ¼ 1, …, Ip; ∀ r

i¼1 H<sup>g</sup><sup>i</sup>

where Nf and Np are the number of unknowns to be determined for neutron flux and delayed neutron precursors, respectively. Substituting expressions (3) and (4) into (1) and (2), and applying the Galerkin process for spatial discretization, as described in [6], the resulting

!; t � �

!; t � �

∈ Ω � ð � 0; T ; (3)

∈ Ω � ð � 0; T ; (4)

∈ Ω � ð � 0; T ; (6)

Ωujukd r!

Ωvlvmd r!

ΩD<sup>g</sup>∇uj∙∇ukd r!

Ων<sup>g</sup><sup>0</sup> <sup>P</sup><sup>g</sup><sup>0</sup>

Ωvjumd r!

<sup>f</sup> ujukd r!

<sup>f</sup> ulvkd r!

: Table 1 contains the expres-

g jk <sup>¼</sup> <sup>Ð</sup>

g0 !g jk <sup>¼</sup> <sup>Ð</sup> Ω P<sup>g</sup><sup>0</sup> !g <sup>s</sup> ujukd r!

gg0 jk <sup>¼</sup> <sup>χ</sup><sup>g</sup><sup>Ð</sup>

gi jm <sup>¼</sup> <sup>λ</sup>iχ<sup>g</sup><sup>Ð</sup>

ig0 lk ¼ β<sup>i</sup> Ð Ων<sup>g</sup><sup>0</sup> <sup>P</sup><sup>g</sup><sup>0</sup> (5)

ð Þt Cið Þt , g ¼ 1, …, G;

!; t � �

tions in terms of base functions defined over Ω, it is possible to write

uk r !� � ϕg

v<sup>m</sup> r !� � Cm

algebraic system of equations can be expressed in a matrix notation as follows:

G

g0 ¼1 Sg0 !<sup>g</sup>ϕ<sup>g</sup><sup>0</sup> ð Þt

X G

and CiðÞ¼ <sup>t</sup> <sup>C</sup><sup>1</sup>

Matrix Type Dimension Elements <sup>M</sup><sup>f</sup> Mass Nf � Nf mf ,jk <sup>¼</sup> <sup>Ð</sup>

<sup>M</sup><sup>p</sup> Mass Np � Np mp,lm <sup>¼</sup> <sup>Ð</sup>

<sup>K</sup><sup>g</sup> Stiffness Nf � Nf <sup>k</sup>

!<sup>g</sup> Mass Nf � Nf s

Fgg<sup>0</sup> Mass Nf � Nf f

Hgi Mass Nf � Np h

Pig<sup>0</sup> Mass Np � Nf p

g0 ¼1 Fgg<sup>0</sup> ð Þ<sup>t</sup> <sup>ϕ</sup><sup>g</sup><sup>0</sup>

� <sup>X</sup> Nf

� <sup>X</sup> Np

m¼1

dt <sup>ϕ</sup><sup>g</sup>ðÞ¼� <sup>t</sup> <sup>K</sup><sup>g</sup>ϕ<sup>g</sup>ð Þ� <sup>t</sup> <sup>X</sup>

<sup>þ</sup> <sup>1</sup> � <sup>β</sup> � �χ<sup>g</sup>

k¼1

As fully explained in [6] and summarized in [1], a simple NFE element is characterized by the fact that for each node, the function unknowns to be determined are the (00) Legendre moment (average) of the unknown function over each face of the node and the (000) Legendre moment over the node volume. Figure 3(a) shows a physical domain Ω graphically represented after generating an xyz mesh. Figure 3(b) shows a cuboid-type node with directions through the faces: (x) Right, Left; (y) Near, Far; (z) Top, Bottom; and C for the average of the function over the node volume. Taking into consideration the general form to build up nodal schemes [7], the moments of a function (at edges and body) over a node like the one shown in Figure 3(b) can be written for the NFE method RTN-0 (Raviart-Thomas-Nédélec).

In the NFE method RTN-0, the normalized zero-order Legendre polynomials defined over the unit cell Ωijk = [�1,+1] � [�1,+1] � [�1,+1] and correlated to each physical cell Ω<sup>e</sup> = Ωijk = [xi ,xi+1] � [yj ,yj+1] � [zk,zk+1] are used to calculate the elements of the matrices in Eqs. (5) and (6).

The matrix elements are quantified introducing the following nodal basis functions [7]:

$$\mathbf{u}\_{\rm L}^{\rm{00}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = -\frac{1}{2} (\mathbf{P}\_{100} - \mathbf{P}\_{200}); \mathbf{u}\_{\rm R}^{\rm{00}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = +\frac{1}{2} (\mathbf{P}\_{100} + \mathbf{P}\_{200});$$

$$\mathbf{u}\_{\rm N}^{\rm{00}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = -\frac{1}{2} (\mathbf{P}\_{010} - \mathbf{P}\_{020}); \mathbf{u}\_{\rm F}^{\rm{00}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = +\frac{1}{2} (\mathbf{P}\_{010} + \mathbf{P}\_{020}); \tag{7}$$

$$\mathbf{u}\_{\rm B}^{\rm{00}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = -\frac{1}{2} (\mathbf{P}\_{001} - \mathbf{P}\_{020}); \mathbf{u}\_{\rm T}^{\rm{00}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = +\frac{1}{2} (\mathbf{P}\_{001} + \mathbf{P}\_{002});$$

$$\mathbf{u}\_{\rm C}^{\rm{000}}(\mathbf{x}, \mathbf{y}, \mathbf{z}) = \mathbf{P}\_{000} - \mathbf{P}\_{200} - \mathbf{P}\_{020} - \mathbf{P}\_{002};$$

where Plpqð Þ¼ x; y; z Plð Þx Ppð Þy Pqð Þz .

Figure 3. Discretization of reactor volume Ω and a local node Ωe. (a) Domain Ω. (b) Physical local node Ωe.

An extensive discussion on nodal diffusion methods can be found in Ref. [7] for space discretization using simplification approaches for calculating the moments over a node.

#### 3.4. Discretization of the time variable

Once the spatial discretization is done, the θ-method can be applied [6] for the discretization of the time variable appearing in the algebraic system given by (5) and (6). For the time integration over the interval (0, T], this interval is divided in L time-steps [tl, tl+1], and the following approach is assumed:

$$\int\_{\mathfrak{h}}^{\mathfrak{h}\_{l+1}} f(\mathbf{t})dt \cong \mathbf{h}\_l \left[\boldsymbol{\Theta}f\_{l+1} + (\mathbf{1} - \boldsymbol{\Theta})f\_l\right] \tag{8}$$

Theoretically, it would be best to determine the flux level which resulted in a critical reactor ð Þ eigenvalue λ<sup>0</sup> ¼ 1 . This could be accomplished by coupling of the NK model with the TH model of the whole reactor. In practice, however, the scaling factor ϕnorm is determined such that the total generated thermal power corresponds to some user-specified value Pth,tot. Before showing how this is done, the relation between the fluxes and the generated thermal power is described. For a given discretization of the xy-plane with pieces of area Δa = Δx�Δy, the thermal

<sup>f</sup> is the volumetric heat generation rate in the fuel in units of [W/cm<sup>3</sup>

differential fuel volume, and the limits zb and zt refer to the coordinates of the bottom and top of the reactor core, respectively. For a given area Δa, the volumetric heat generation rate q<sup>000</sup>

G

g0 ¼1 Σg0 <sup>f</sup> ð Þ<sup>z</sup> <sup>ϕ</sup><sup>g</sup><sup>0</sup>

where ϕnorm is a dimensionless factor, Efiss is the energy released by a nuclear fission reaction in

ðzt

X G

g0 ¼1 Σg0 <sup>f</sup> ð Þ<sup>z</sup> <sup>ϕ</sup><sup>g</sup><sup>0</sup>

zb

In a more general way, for a reactor volume V composed by the union of sub-volumes Ve (see

Ne

X G

> 3 5 �1

<sup>c</sup> <sup>¼</sup> <sup>ϕ</sup>normϕ<sup>e</sup>

g0 ¼1 Σg0 <sup>f</sup> ,eϕ<sup>g</sup><sup>0</sup>

e¼1

X G

g0 ¼1 κg0 <sup>f</sup> ,eϕ<sup>g</sup><sup>0</sup> <sup>e</sup> Ve

Therefore, using the reference total thermal power specified by the code user, the flux normal-

Ne

2 4

e¼1

<sup>f</sup> ,e <sup>¼</sup> EfissΣ<sup>g</sup><sup>0</sup>

calculated as above, the actual thermal power distributions in the reactor core can be calculated

introduce the value of Efiss. This value is used as an average energy released of �200 MeV

<sup>f</sup> ð Þz da � dz, dV ¼ da � dz; (10)

ð Þz ; (11)

ð Þz da � dz: (12)

<sup>e</sup> Ve: (13)

; (14)

. Nevertheless, it is necessary to

<sup>f</sup> ,e: With the flux normalization factor ϕnorm

], dV is a

Nuclear Reactor Simulation

13

http://dx.doi.org/10.5772/intechopen.79723

∙s)]. Thus,

<sup>f</sup> ð Þz

power Pth,tot can be expressed as follows:

where q<sup>000</sup>

Eq. (10) is written as

ization factor can be written as

where the factors "kappa-fission" are κ<sup>g</sup><sup>0</sup>

using the current neutron flux in the reactor core ϕ<sup>e</sup>

Pth,tot <sup>¼</sup> <sup>X</sup>

in an elevation z may be written in terms of the fluxes as

Figure 3), the total thermal power can be expressed as

q<sup>000</sup>

Pth,tot <sup>¼</sup> <sup>ϕ</sup>normEfiss<sup>X</sup>

Δa

ð zt

zb q<sup>000</sup>

<sup>f</sup> ð Þ¼ <sup>z</sup> <sup>ϕ</sup>normEfiss <sup>X</sup>

[MeV/fission], and the sum over g<sup>0</sup> is the volumetric fission rate in [fissions/(cm<sup>3</sup>

Δa

Pth,tot <sup>¼</sup> <sup>ϕ</sup>normEfiss<sup>X</sup>

<sup>ϕ</sup>norm <sup>¼</sup> Pth,tot <sup>X</sup>

(i.e.,), based on the energies released by the fission of the U<sup>235</sup> nuclei [8].

where hl <sup>¼</sup> tlþ<sup>1</sup> � tl, f <sup>l</sup> <sup>¼</sup> f tð Þ<sup>l</sup> , f <sup>l</sup>þ<sup>1</sup> <sup>¼</sup> f tð Þ <sup>l</sup>þ<sup>1</sup> , and <sup>θ</sup> is the time integration parameter.

For time integration, parameters θ<sup>f</sup> and θ<sup>p</sup> for neutron flux and delayed neutron precursors are considered with values in the interval [0, 1], giving different time integration schemes [6].

Once the formulation to be used for time integration is established, the NfG + NpI system of equations that was spatially discretized, Eqs. (5) and (6) are discretized over the interval (0,T]. Integrating the referred equations over the time interval [tl, tl+1] using approximation (8), the following set of equations is generated:

$$\mathbf{A}\_{l+1}\boldsymbol{\Phi}\_{l+1} = \mathbf{Q}\_{l};\ l = 0, 1, 2, \dots;$$

$$\boldsymbol{\Phi}\_{l+1} = \begin{bmatrix} \phi\_{l+1}^{\mathrm{I}}, \dots, \phi\_{l+1}^{\mathrm{G}} \end{bmatrix}^{\mathrm{T}};\ \ Q\_{l} = \begin{bmatrix} \mathbf{S}\_{f,l}^{\mathrm{I}}, \dots, \mathbf{S}\_{f,l}^{\mathrm{G}} \end{bmatrix}^{\mathrm{T}};$$

For a known vector Φ<sup>l</sup> the algebraic system (9) is solved for the neutron fluxes Φ<sup>l</sup>þ1. Therefore, the computing process requires an initial flux vector for the first time step, which is used in (9) to determine new neutron fluxes at the end of the time step, thus using these neutron fluxes to calculate a new delayed neutron precursor concentration vector. This process is sequentially performed for each time step over the total time interval (0,T].

### 4. Reactor power distribution

Once the computer model to solve the reactor kinetics Eqs. (1) and (2) is able to provide the neutron flux profile, the next objective is to know the power distribution in the reactor configuration. It is necessary to be aware that the neutron flux is by itself the shape of the power distribution in multiplicative materials. The numerical methods presented in previous sections to solve Eq. (9) produce an algorithm capable to obtain the neutron flux profile for a reactor steady state. The calculated neutron flux has the following property over the domain Ω: ϕ � � � � ¼ 1. To determine the real average neutron flux in the reactor core, ϕc, it is necessary to specify the magnitude of the fluxes. For instance, a flux normalization factor ϕnorm can be introduced such that <sup>ϕ</sup><sup>c</sup> <sup>¼</sup> <sup>ϕ</sup>norm<sup>ϕ</sup> neutrons cm2∙seg h i:

Theoretically, it would be best to determine the flux level which resulted in a critical reactor ð Þ eigenvalue λ<sup>0</sup> ¼ 1 . This could be accomplished by coupling of the NK model with the TH model of the whole reactor. In practice, however, the scaling factor ϕnorm is determined such that the total generated thermal power corresponds to some user-specified value Pth,tot. Before showing how this is done, the relation between the fluxes and the generated thermal power is described. For a given discretization of the xy-plane with pieces of area Δa = Δx�Δy, the thermal power Pth,tot can be expressed as follows:

An extensive discussion on nodal diffusion methods can be found in Ref. [7] for space discretization using simplification approaches for calculating the moments over a node.

Once the spatial discretization is done, the θ-method can be applied [6] for the discretization of the time variable appearing in the algebraic system given by (5) and (6). For the time integration over the interval (0, T], this interval is divided in L time-steps [tl, tl+1], and the following

<sup>f</sup>ð Þ<sup>t</sup> dt ffi hl <sup>θ</sup><sup>f</sup> <sup>l</sup>þ<sup>1</sup> <sup>þ</sup> ð Þ <sup>1</sup> � <sup>θ</sup> <sup>f</sup> <sup>l</sup>

For time integration, parameters θ<sup>f</sup> and θ<sup>p</sup> for neutron flux and delayed neutron precursors are considered with values in the interval [0, 1], giving different time integration schemes [6].

Once the formulation to be used for time integration is established, the NfG + NpI system of equations that was spatially discretized, Eqs. (5) and (6) are discretized over the interval (0,T]. Integrating the referred equations over the time interval [tl, tl+1] using approximation (8), the

A<sup>l</sup>þ<sup>1</sup>Φ<sup>l</sup>þ<sup>1</sup> ¼ Ql; l ¼ 0, 1, 2, …, ; (9)

f ,l ; …; SG f ,l h i<sup>T</sup>

; Ql <sup>¼</sup> <sup>S</sup><sup>1</sup>

where hl <sup>¼</sup> tlþ<sup>1</sup> � tl, f <sup>l</sup> <sup>¼</sup> f tð Þ<sup>l</sup> , f <sup>l</sup>þ<sup>1</sup> <sup>¼</sup> f tð Þ <sup>l</sup>þ<sup>1</sup> , and <sup>θ</sup> is the time integration parameter.

<sup>l</sup>þ<sup>1</sup>; …; <sup>ϕ</sup><sup>G</sup>

� �<sup>T</sup>

lþ1

For a known vector Φ<sup>l</sup> the algebraic system (9) is solved for the neutron fluxes Φ<sup>l</sup>þ1. Therefore, the computing process requires an initial flux vector for the first time step, which is used in (9) to determine new neutron fluxes at the end of the time step, thus using these neutron fluxes to calculate a new delayed neutron precursor concentration vector. This process is sequentially

Once the computer model to solve the reactor kinetics Eqs. (1) and (2) is able to provide the neutron flux profile, the next objective is to know the power distribution in the reactor configuration. It is necessary to be aware that the neutron flux is by itself the shape of the power distribution in multiplicative materials. The numerical methods presented in previous sections to solve Eq. (9) produce an algorithm capable to obtain the neutron flux profile for a reactor steady state. The calculated neutron flux has the following property over the domain Ω:

� ¼ 1. To determine the real average neutron flux in the reactor core, ϕc, it is necessary to specify the magnitude of the fluxes. For instance, a flux normalization factor ϕnorm can be

> cm2∙seg h i :

� � (8)

;

ðtlþ<sup>1</sup> tl

<sup>Φ</sup><sup>l</sup>þ<sup>1</sup> <sup>¼</sup> <sup>ϕ</sup><sup>1</sup>

performed for each time step over the total time interval (0,T].

3.4. Discretization of the time variable

following set of equations is generated:

4. Reactor power distribution

introduced such that <sup>ϕ</sup><sup>c</sup> <sup>¼</sup> <sup>ϕ</sup>norm<sup>ϕ</sup> neutrons

ϕ � � �

approach is assumed:

12 New Trends in Nuclear Science

$$\mathbf{P\_{th,tot}} = \sum\_{\Delta a} \int\_{\mathbf{z}\_b}^{\mathbf{z}\_l} \mathbf{q\_f'''(z)} \mathbf{da} \cdot \mathbf{dz}, \text{ dV} = \mathbf{da} \cdot \mathbf{dz} \tag{10}$$

where q<sup>000</sup> <sup>f</sup> is the volumetric heat generation rate in the fuel in units of [W/cm<sup>3</sup> ], dV is a differential fuel volume, and the limits zb and zt refer to the coordinates of the bottom and top of the reactor core, respectively. For a given area Δa, the volumetric heat generation rate q<sup>000</sup> <sup>f</sup> ð Þz in an elevation z may be written in terms of the fluxes as

$$q\_f'''(z) = \phi\_{norm} E\_{\text{fits}} \sum\_{\mathbf{g}'=1}^G \Sigma\_f^{\mathbf{g}'}(z) \phi^{\mathbf{g}'}(z);\tag{11}$$

where ϕnorm is a dimensionless factor, Efiss is the energy released by a nuclear fission reaction in [MeV/fission], and the sum over g<sup>0</sup> is the volumetric fission rate in [fissions/(cm<sup>3</sup> ∙s)]. Thus, Eq. (10) is written as

$$P\_{th,tot} = \phi\_{nmm} E\_{fks} \sum\_{\Delta a} \int\_{z\_b}^{z\_l} \sum\_{\mathbf{g}'=1}^{G} \Sigma\_f^{\mathbf{g}'}(z) \phi^{\mathbf{g}'}(z) da \cdot dz. \tag{12}$$

In a more general way, for a reactor volume V composed by the union of sub-volumes Ve (see Figure 3), the total thermal power can be expressed as

$$P\_{th,tot} = \phi\_{norm} E\_{f\text{fss}} \sum\_{\varepsilon=1}^{N\_{\varepsilon}} \sum\_{\mathbf{g}'=1}^{G} \Sigma\_{f,\varepsilon}^{\mathbf{g}'} \phi\_{\varepsilon}^{\mathbf{g}'} V\_{\varepsilon}.\tag{13}$$

Therefore, using the reference total thermal power specified by the code user, the flux normalization factor can be written as

$$\phi\_{norm} = P\_{th,tot} \left[ \sum\_{e=1}^{N\_e} \sum\_{\mathbf{g}'=1}^G \mathbf{k}\_{f,e}^{\mathbf{g}'} \phi\_e^{\mathbf{g}'} V\_e \right]^{-1};\tag{14}$$

where the factors "kappa-fission" are κ<sup>g</sup><sup>0</sup> <sup>f</sup> ,e <sup>¼</sup> EfissΣ<sup>g</sup><sup>0</sup> <sup>f</sup> ,e: With the flux normalization factor ϕnorm calculated as above, the actual thermal power distributions in the reactor core can be calculated using the current neutron flux in the reactor core ϕ<sup>e</sup> <sup>c</sup> <sup>¼</sup> <sup>ϕ</sup>normϕ<sup>e</sup> . Nevertheless, it is necessary to introduce the value of Efiss. This value is used as an average energy released of �200 MeV (i.e.,), based on the energies released by the fission of the U<sup>235</sup> nuclei [8].

In summary, once the NK model is used to generate the neutron flux distribution in the reactor core, expression (12) can be used to calculate the thermal power being generated along all the nodes in a thermal-hydraulic channel of area Δa and height H. This thermal power can be the axial power profile needed by the TH model to produce the thermal-hydraulic state corresponding to the generated thermal power.

condition for each successive ΔT. Achieving converge for each ΔT with respective reactor core conditions means to produce a time-dependent behavior of the reactor condition over the total

The TH model comprises the solution of the mass, momentum, and energy conservation equations in the three regions contemplated by the channel: (1) one phase, (2) subcooled boiling, and (3) bulk boiling. The system receives heat through a non-uniform source whose profile is axially defined plane by plane. This axial use of the power profile allows the inclusion of a wide range of axial profiles, from relatively flat to profiles with their peak value at some

In the following subsections, there are several expressions for which the corresponding param-

The heat transfer and temperature distribution in the fuel and cladding can be calculated by a simple model where the heat diffusion equation is solved in one dimension (radial) for a fuel rod, since the conduction in axial direction is small compared to the radial one, it can be

> ð Þ� t 1 R0 g

Tf ð Þ� <sup>t</sup> Tc ð Þ<sup>t</sup> � <sup>1</sup>

to the refrigerant fluid is calculated by the Dittus-Boelter or Chen correlation, depending on the type of flow, which can be in one or two phases. These equations are used for the radial

The conservation equations of mass, energy, and momentum are applied in this case to a flow of water along a vertical channel, where the dynamics of the fluid heated by the wall of the fuel

∂Gm

<sup>∂</sup><sup>z</sup> � <sup>f</sup>Gmj j Gm 2Der<sup>m</sup>

> ∂p ∂z þ

∂r<sup>m</sup> ∂t þ

¼ � <sup>∂</sup><sup>p</sup>

<sup>00</sup>Ph Az þ ∂p ∂t þ Gm rm

R0 c

<sup>c</sup> represent thermal resistances per unit length. The coefficient of heat transfer

Tfð Þ� <sup>t</sup> Tcð Þ<sup>t</sup> (15)

<sup>∂</sup><sup>z</sup> <sup>¼</sup> <sup>0</sup> (17)

fGmj j Gm 2Der<sup>m</sup>

� rmgcosθ (18)

(19)

Tc ð Þ� <sup>t</sup> Tm ð Þ<sup>t</sup> (16)

Nuclear Reactor Simulation

15

http://dx.doi.org/10.5772/intechopen.79723

time interval T.

axial point in each channel in the core.

neglected. An energy balance per unit length yields

mccpc

averaging of the temperatures in the fuel rod.

mf cpf

dTc dt <sup>¼</sup> <sup>1</sup> R0 g

is modeled. Conservation equations can be expressed as [10]

∂ ∂z

G2 m r<sup>þ</sup> m 

∂hm <sup>∂</sup><sup>z</sup> <sup>¼</sup> <sup>q</sup>

∂Gm ∂t þ

rm ∂hm <sup>∂</sup><sup>t</sup> <sup>þ</sup> Gm dTf dt <sup>¼</sup> <sup>q</sup><sup>0</sup>

eters are defined in Refs. [10, 11].

5.1. Heat transfer in the fuel

where R<sup>0</sup>

<sup>g</sup> and R<sup>0</sup>

5.2. Reactor coolant dynamics
