Effective point kinetic parameters calculation in Tehran research reactor using deterministic and probabilistic methods

NUCLEAR ENERGY SCIENCE AND ENGINEERING

Effective point kinetic parameters calculation in Tehran research reactor using deterministic and probabilistic methods

M. Kheradmand Saadi，

A. Abbaspour

Nuclear Science and Techniques

Vol.28, No.12

Article number 171

Published in print 01 Dec 2017

Available online 16 Nov 2017

DOI：10.1007/s41365-017-0323-7

38601

The exact calculation of point kinetic parameters is very important in nuclear reactor safety assessment and most sophisticated safety codes such as RELAP5, PARCS, DYN3D, PARET, are using these parameters in their dynamic models. These parameters include effective delayed neutron fractions as well as mean generation time. These parameters are adjoint weighted, and adoint flux is employed as a weighting function in their evaluation. Adjoint flux calculation is an easy task for most of deterministic codes, but its evaluation is cumbersome for Monte Carlo codes. However, in recent years, some sophisticated techniques have been proposed for Monte Carlo based point kinetic parameters calculation without any need of adjoint flux. The most straightforward scheme is known as the "Prompt Method" and has been used widely in literature. The main objective of this article is dedicated to point kinetic parameters calculation in Tehran Research Reactor (TRR) using deterministic as well as probabilistic techniques. WIMS-D5B and CITATION codes have been used in deterministic calculation of forward and adjoint fluxes in the TRR core. On the other hand, the MCNP Monte Carlo code has been employed in the "Prompt Method" scheme for effective delayed neutron fraction evaluation. Deterministic results have been cross checked with probabilistic ones and validated with SAR and experimental data. In comparison with experimental results, the relative differences of deterministic as well as probabilistic methods are 7.6% and 3.2%, respectively. These quantities are 10.7% and 6.4% respectively in comparison with SAR report.

Point kinetic parametersTehran research reactorAdjoint fluxPrompt methodDeterministic methodProbabilistic method

1 Introduction

One of the usual simplifications, which is used widely in nuclear reactor analysis, is the neutron diffusion equation. It considers a linear angular dependency of neutron flux and scattering cross section in the neutron transport equation. The numerical solution of diffusion equation, especially in time dependent cases, is very difficult yet. However, some special efforts have been developed to simplify the time dependent neutron diffusion equation. The most famous treatment, which focused on space variable elimination, is known as the “Point Kinetic Model”. This model assumes that the spatial dependency of neutron flux can be described by a single spatial mode (fundamental mode). Moreover, this shape function is time independent and will be constant during each transient ^[1]. By this assumption, the time dependent neutron flux is divided into two independent functions: a time dependent amplitude function and a time independent shape function ^[1]_,

φ (r, t) = ϑ n (t) ψ (r),

(1)

where, $ϑ$ is neutron velocity, n is neutron density and $ψ$ is shape function and dimensionless in this formula. This technique can be employed in one group time dependent neutron diffusion equation, including delayed neutrons and point kinetic equations, and can be derived as follows ^[1]:

\begin{array}{l} \frac{d n}{d t} = [\frac{ρ (t) - β}{Λ}] n (t) - \sum_{i = 1}^{6} λ_{i} C_{i}, \\ \frac{d C_{i}}{d t} = \frac{β_{i}}{Λ} n (t) - λ_{i} C_{i} i = 1, \dots, 6 . \end{array}

(2)

These equations are a set of seven coupled ordinary differential equations in time and describe both the time dependency of the neutron population in a reactor and the decay of delayed neutron precursors. As usual, delayed neutron precursors are presented in six groups. In these equations, $n (t)$ is neutron population density; $C_{i} (t)$ is nucleus precursor density of i^th group; β is the effective delayed neutron fraction and usually presented as β_eff; Λ is the mean generation time; λ_i is the decay constant of i^th group, and β_i is delayed neutron fraction of i^th group. From these parameters, β_i, β_eff, and Λ are spectrum and adjoint weighted and may vary from reactor to reactor based on different compositions and spectrums. These parameters - which are known as point kinetic parameters - are very important in the analytical safety analysis of each reactor. The reactivity is directly related to the delayed neutron fraction and determines the prompt criticality. Consequently, an exact calculation of point kinetic parameters is very important in each nuclear reactor design. Most of sophisticated safety analysis codes such as RELAP5, PARCS, DYN3D, and PARET are using these parameters in their dynamic models.

Traditionally, the deterministic approach has been employed to evaluate the effective point kinetic parameters. This involve an adjoint and spectrum weighting of the delayed neutron production rate and hence requires a calculation of both forward as well as adjoint fluxes and, on top of that, a suitable post-processing to calculate the weighted production rate ^[2]. Due to the free availability and versatility of WIMSD and CITATION codes, they have been used widely by developing countries for deterministic calculation of point kinetic parameters. Some efforts have been performed to evaluate the point kinetic parameters in PARR-1 as a 10 MW pool type research reactor ^[3-5]. In these calculations, forward as well as adjoint fluxes were evaluated based on the microscopic cross section approach. WIMSD/4 calculates microscopic cross sections and an auxiliary program, known as BORGES, writes them in the format of the CITATION input requirement. Calculations are in very good agreement with experimental and PARR-1 FSAR data.

On the other hand, adjoint calculations are cumbersome in continuous energy Monte Carlo codes ^[6]. However, in recent years, some techniques have been proposed to evaluate the Monte Carlo based point kinetic parameters by formulations which do not require the adjoint flux. The simplest and first order method in the Monte Carlo evaluation of point kinetic parameters is known as the “Prompt Method” scheme, which has been described by Meulekamp and Van der Marck ^[7]. The results of this method are in good agreement with our experiment and other adjoint weighted methods and usually is known as an acceptable method in Monte Carlo based delayed neutron calculation ^[8].

The main objective of this article is dedicated to the calculation of point kinetic parameters in the Tehran Research Reactor (TRR). Deterministic, as well as probabilistic approaches, have been employed and compared with each other. The use of macroscopic cross sections in deterministic calculations and the Monte Carlo "Prompt Method" technique in the probabilistic approach are main innovative aspects in this work. The need of adjoint flux as a weighting function and the meaning of "effective" will be explained in Sect.2. The Tehran Research Reactor by considering different fuel element types will be described in Sect. 3. The deterministic, as well as probabilistic, methodologies in effective point kinetic parameters evaluation will be described in Sect. 4. Results and discussion will be presented in Sect. 5 and finally, the article will be concluded in Sect. 6.

2 Effective delayed neutron parameters

Essentially, most fission neutrons appear instantaneously after the fission event. These neutrons are referred to as prompt neutrons. However, a very few fraction of neutrons are appeared with an appreciable delayed time from the subsequent decay of radioactive fission products. Although only a few fractions of fission neutrons are delayed, they are vital for the effective control of the fission chain reaction. Delayed neutrons do not have the same properties as prompt ones released directly from fission. The averaged energy of prompt neutrons is about 2 MeV, which is much greater than the averaged energy of delayed neutrons (about 0.5 MeV). This fact that delayed neutrons are born at lower energies has two significant impacts on the way they precede through the neutron life cycle. Firstly, delayed neutrons have a much lower probability to cause fast fissions in comparison with prompt ones because their averaged energy is less than the threshold required for fast fission. Secondly, delayed neutrons have a lower probability to leak out of the core because they are born at lower energies and subsequently travel shorter distances in comparison with fast neutrons. In other words, the delayed and prompt neutrons have a difference in their effectiveness in producing a subsequent fission event. Since the energy distribution of delayed neutrons is different from group to group, different groups of delayed neutrons have a different effectiveness.

To deal with this situation, it is necessary to define an important function, $ϕ^{†} (r, E)$ , which is the probability that a neutron introduced at position r and energy E will ultimately results in fission. Then, the relative importance (to the production of a subsequent fission) of delayed and prompt neutrons in group i of isotope q are $I_{d i}^{q}$ and $I_{p}^{q}$ , respectively ^[9].

I_{d i}^{q} = \int d V \int_{o}^{\infty} d E χ_{d i}^{q} (E) ϕ^{†} (r, E) \int_{0}^{\infty} d E^{'} ν σ_{f}^{q} (E^{'}) N_{q} (r) ϕ (r, E^{'})

(3)

I_{p}^{q} = \int d V \int_{o}^{\infty} d E χ_{p}^{q} (E) ϕ^{†} (r, E) \int_{0}^{\infty} d E^{'} ν σ_{f}^{q} (E^{'}) N_{q} (r) ϕ (r, E^{'})

(4)

The main variables in these equations are as follows: $χ_{di}^{q}$ : Delayed neutron spectrum of the fissionable isotope, q, in group i $χ_{p}^{q}$ : Prompt neutron spectrum of the fissionable isotope, q $ϕ^{†}$ : Adjoint flux distribution $ν$ : Neutron fission yield $σ_{f}^{q}$ : Microscopic fission cross section of the fissionable isotope, q $N_{q}$ : Atom density of the fissionable isotope, q $I_{di}^{q}$ : Relative importance of delayed neutrons in group i of isotope q $I_{p}^{q}$ : Relative importance of prompt neutrons of isotope q

For the fissionable isotope, q, the relative effective delayed neutrons yield in the i^th group of the delayed neutrons is $I_{di}^{q} β_{i}^{q}$ . $β_{i}^{q}$ is the group i of delayed neutrons yield for the fissionable isotope, q. Consequently the effective delayed neutrons fraction of isotope q in group i is $\frac{I_{di}^{q} β_{i}^{q}}{I_{P}^{q}}$ ^[9]:

β_{i}^{q - e f f} = \frac{\int d V \int_{0}^{\infty} d E ϕ^{†} (r, E) χ_{d i}^{q} (E) β_{i}^{q} \int_{0}^{\infty} υ Σ_{f}^{q} (r, E^{'}) ϕ (r, E^{'}) d E^{'}}{\int d V \int_{0}^{\infty} d E ϕ^{†} (r, E) χ_{p}^{q} (E) \int_{0}^{\infty} υ Σ_{f}^{q} (r, E^{"}) Φ (r, E^{'}) d E^{'}} .

(5)

The discretized format of Eq. (5) is as follow:

β_{i}^{q - e f f} = \frac{\sum_{n = 1}^{N_{m e s h}} V_{n} \sum_{g = 1}^{G} ϕ_{g, n}^{†} χ_{d i}^{g} β_{i}^{q} \sum_{g^{'} = 1}^{G} υ Σ_{f, g^{'}, n} Φ_{g^{'}, n}}{\sum_{n = 1}^{N_{m e s h}} V_{n} \sum_{g = 1}^{G} ϕ_{g, n}^{†} χ_{p}^{g} \sum_{g^{'} = 1}^{G} υ Σ_{f, g^{'}, n} Φ_{g^{'}, n}} .

(6)

G is the total number of groups and V_n is the volume of the n^th mesh point. In a mixture of fissionable isotopes, the effective delayed neutron fraction in the i^th group ( $β_{i}^{eff}$ ) becomes the average of different $β_{i}^{q - eff}$ weighted by amount of fission due to each isotope ^[10]. The total effective delayed neutron fraction will be the sum of these fractions at different groups:

β_{e f f} = \sum_{i = 1}^{6} β_{i}^{e f f} .

(7)

Such a manner should be repeated for mean generation time calculation, and finally, the following formula will be obtained ^[10]:

Λ = \frac{\int d V \int_{0}^{\infty} d E ϕ^{†} (r, E) [\frac{1}{v (E)}] ϕ (r, E)}{\int d V \int_{0}^{\infty} d E ϕ^{†} (r, E) χ_{p} (E) \int_{0}^{\infty} υ Σ_{f} (r, E^{"}) Φ (r, E^{'}) d E^{'}}

(8)

The discretized format of Eq. (8) is as follow:

Λ = \frac{\sum_{n = 1}^{N_{m e s h}} V_{n} \sum_{g = 1}^{G} ϕ_{g, n}^{†} [\frac{1}{v^{g}}] ϕ_{g, n}}{\sum_{n = 1}^{N_{m e s h}} V_{n} \sum_{g = 1}^{G} ϕ_{g, n}^{†} χ_{p}^{g} \sum_{g^{'} = 1}^{G} υ Σ_{f, g^{'}, n} Φ_{g^{'}, n}} \cdot

(9)

3 Tehran research reactor

The Tehran Research Reactor (TRR) is a 5 megawatt pool-type light-water moderated, heterogeneous solid fuel reactor in which the water is also used for cooling and shielding ^[11]. The TRR core is immersed in either section of a two-section, concrete pool filled with water. The utilization of the reactor is essentially for research, training, and production of radioisotopes ^[11]. The reactor core is composed of MTR-type fuel assemblies inserted in the grid plate. The assemblies may be arranged in a variety of lattice patterns depending on the experimental requirements. The main characteristics of the TRR core have been shown in Table 1. A series of fifty-four holes, capable of accommodating the end fittings of fuels, are arranged in a 9 × 6 rectangular lattice. The first core configuration is a 19 fuel element arrangement reflected by water, as shown in Fig.1. The original core of TRR was High Enriched Uranium (HEU) fuels. Upon procurement of the next cores, Low Enriched Uranium (LEU) fuels had to be considered ^[11]. There are two types of LEU fuel elements named as the Standard Fuel Element (SFE) and Control Fuel Element (CFE). Each fuel type and the core physical characteristics are described in following subsections.

Main characteristics of Tehran Research Reactor ^[11]

Fuel elements:
U²³⁵ per Standard Fuel Element (SFE)	214 gr
U²³⁵ per Control Fuel Element (CFE)	76 gr
U per fuel plate	290 gr
Number of plates per fuel element	19 for SFE
	14 for CFE
Meat:
Enriched U₃O₈	20% in weight of U²³⁵
U density	2.9617 gr/cm³
Meat density	4.76 gr/cm³
Void fraction	10.0%
Weight percentage	U²³⁵12.45%, U²³⁸ 49.78%, O 11.18%, Al 26.59%
Aluminum meat	Purity 99.6%
	Density=2.7 gr/cm³
Shim and safety rods absorber	Four Ag-In-Cd Alloy (80%,15%,5% in weight respectively)Density 10.17 gr/cm³
Control rods’ cladding material	AISI-316/L stainless steel
	Density= 7.95 gr/cm³
Regulating rod	One AISI-316/L stainless steelDensity= 7.95 gr/cm³
Grid plate	Grid array X-Y
	Pitch: 7.71×8.1 cm
Coolant/Moderator	Light water/ Light water
Reflectors	Light water/Graphite

Fig.1.

(Color online) TRR First Core Configuration

3.1 Standard fuel element (SFE)

The SFE is 20% enriched in weight of ²³⁵U and has 19 flat fuel plates inserted in two grooved side plates (lateral walls), as shown in Fig.2 and Fig.3. The meat is made of U₃O₈ powder dispersed in a pure aluminum matrix. Each fuel plate is made up of fuel meat and cladding which seals it off hermetically while isolating it from the coolant. Cladding consists of a frame and two covers in an annealed aluminum Al-6061 alloy. The gap between the fuel plates and vertical through hole in the end fitting assures proper cooling water flow during the operation. Side plates and external plates are fastened to the end fitting at the assembly's lower end by means of TIG welding. The main and geometric data of SFE is summarized in Table 2.

Summary of Main and Geometric Data for LEU-SFE ^[11]

Enrichment	20 %
Number of fuel plates	19
SFE dimensions	8.01×7.6×91.8 cm
Plate thickness	0.15 cm
Clad thickness	0.04 cm
Water channel thickness	0.27 cm
Meat thickness	0.07 cm
Meat width	6.0 cm
Meat length	61.5
Lateral wall thickness (side plate)	0.45 cm
Inlet/exit channel entrance length	4.55 cm
Inner diameter of inlet/exit Nozzle	5.30 cm
Outer diameter of inlet/exit Nozzle	6.16 cm
Coolant flow area	33.92 cm²
Heat transfer area	14022.0 cm²
Meat material	U₃O₈-Al
Fuel plate cladding and side wall material	Al-6061

Fig.2.

(Color online) Standard Fuel Element (unit: mm)

Fig.3.

(Color online) Standard Fuel Element in 1/4 Symmetry (unit: mm)

3.2 Control fuel element (CFE)

The control fuel element (CFE) has 14 flat fuel plates and accommodates fork-type control rods in a lateral position of fuel assembly as shown in Fig. 4 and Fig. 5. Positioning four shim safety rods and one regulating rod into the reactor core controls the reactor. The shim safety control rod is composed of two reactivity control plates (absorbing plates) and knuckle sub-assembly. The absorber plate is an alloy of Silver, Indium, and Cadmium (80%, 15%, and 5% wt. respectively), while the regulating rod is made of stainless steel. All absorbing rods are fork type. Shim rods facilitate, the start-up and operation of the TRR and ensure safe shut down of the reactor at any moment, while the regulating rod is used for fine reactivity insertion either automatically or manually. Control rods are dropped into the core upon receiving a trip signal from various safety channels. The main and geometrical data of the CFE is summarized in Table 3.

Summary of Main and Geometric Data for CFE ^[11]

Enrichment	20 %
Number of fuel plates	14
CFE dimension	8.01×7.6×161.4 cm
Plate thickness	0.15 cm
Clad thickness	0.04 cm
Water channel thickness	0.27 cm
Meat thickness	0.07 cm
Meat width	6.0 cm
Meat length	61.5
Lateral wall thickness (side plate)	0.45 cm
Inlet/exit channel entrance length	4.55 cm
Inner diameter of inlet/exit Nozzle	5.30 cm
Outer diameter of inlet/exit Nozzle	6.16 cm
Coolant flow area	25.81 cm²
Heat transfer area	10332 cm²
Fuel plate cladding and side walls material	Al-6061
Meat material	U₃O₈-Al
Fuel plate cladding and side walls material	Al-6061
Absorber material for shim safety rods	Ag-In-Cd
Absorber material for fine regulating rod	AISI-316L SS

Fig.4.

Control Fuel Element (unit: mm)

Fig.5.

(Color online) Control Fuel Element in 1/4 Symmetry (unit: mm)

4 Calculation methodologies

4.1 Deterministic approach

Deterministic calculations of point kinetic parameters in the TRR have been performed using WIMS-D5B ^[12] and CITATION-LDI2 ^[13] codes for cell and core calculations, respectively. In order to distinguish between prompt and delayed neutron spectrums, the total number of energy groups in these calculations should be more than the usual two groups' structure. Since in later one fractions of neutrons, which are born in fast and thermal groups, are one and zero, respectively. In other words, in the two groups calculation, all of the neutrons (either prompt or delayed) are born in fast group, and hence, it's not possible to distinguish between prompt and delayed neutrons. The six groups' structure has been considered for cell calculation. The bounds of each energy group - based on the 69 groups' structure of WIMS library - as well as the prompt fission spectrum are shown in Table 4. First two groups are based on the fast group structure of the WIMS library. The first group is above the threshold fission of ²³⁸U and the second one is considered as a reminder. Third and fourth groups are considered for the resonance region of the library. And finally, the last two groups are based on the thermal group structure.

The energy bounds in six group energy discretization

Energy Group	Energy boundariesbased on WIMS library	Prompt fission spectrum
1	10.000 Mev - 2.231 Mev	3.564777E-01
2	2.231 Mev - 9.118 kev	6.431600E-01
3	9.118 kev - 906.898 ev	3.488916E-04
4	906.00 ev - 4.0000 ev	1.339680E-05
5	4.0000 ev - 0.0250 ev	0.000000E+00
6	0.02500 ev - 0.0000 ev	0.000000E+00

If one assumes a Maxwellian distribution, the spectrum of delayed neutrons ( $χ_{d}$ ) is calculated by the following expression ^[14]:

χ_{d i}^{g} = \int_{G_{1}}^{G_{2}} 2.07296489 {(\bar{E})}^{- 1.5} \sqrt{E} e^{(- \frac{1.5 E}{\bar{E}})} d E,

(10)

where

$χ_{d i}^{g}$

is the fraction of the i^th group of delayed neutrons in the g^th energy group,

$\bar{E}$

is the Mean energy of the i^th group of delayed neutrons, and G1, G2 are the energy boundaries of the g^th energy group.

The delayed neutron spectrum of each group has been calculated and shown in Table 5. The mean energy of each delayed neutron group has been taken from Ref.[14]. On the other hand, bounding energies are based on six groups' structure in cell calculation and taken from Table 4. Also, delayed neutron parameters of ²³⁵U and ²³⁸U have been presented in Table 6 and Table 7, respectively.

Energy and delayed neutron spectrum of each delayed group

Family number	$\bar{E}$ (MeV) ^[14]	Energy group	Spectrum
1	0.25	1	6.5722e-06
		2	9.9070e-01
		3	9.0000e-03
		4	3.0087e-04
		5	8.8402e-08
		6	4.3702e-11
2	0.46	1	7.6000e-03
		2	9.8960e-01
		3	2.7000e-03
		4	8.9908e-05
		5	2.6369e-08
		6	1.3036e-11
3	0.405	1	1.4000e-03
		2	9.9440e-01
		3	4.1000e-03
		4	1.3356e-04
		5	3.9189e-08
		6	1.9373e-11
4	0.45	1	1.2900e-02
		2	9.8470e-01
		3	2.4000e-03
		4	7.7188e-05
		5	2.2635e-08
		6	1.1190e-11
5	0.42	1	1.2000e-03
		2	9.9450e-01
		3	4.2000e-03
		4	1.3835e-04
		5	4.0597e-08
		6	2.0070e-11
6	0.18	1	4.2082e-08
		2	9.8490e-01
		3	1.4600e-02
		4	4.9185e-04
		5	1.4470e-07
		6	7.1532e-11

Delayed neutron parameters of ²³⁵U ^[9]

Delayed neutron group	Fast Neutrons		Thermal Neutrons
	Decay constantλ_i (S^-1)	Delayed neutron yield (β=0.0064)	Decay constantλ_i (S^-1)	Delayed neutron yield (β=0.0067)
1	0.0127	0.0002432	0.0124	0.0002211
2	0.0317	0.0013632	0.0305	0.0014673
3	0.115	0.0012032	0.111	0.0013132
4	0.311	0.0026048	0.301	0.0026465
5	1.40	0.0008192	1.14	0.0007705
6	3.87	0.0001664	3.01	0.0002814

Delayed neutron parameters of ²³⁸U ^[9]

Delayed Neutron Group	Fast Neutrons
	Decay Constantλ_i (S^-1)	Delayed Neutron Yield (β =0.0164)
1	0.0132	0.0002132
2	0.0321	0.0022468
3	0.139	0.0026568
4	0.358	0.0063632
5	1.41	0.00369
6	4.02	0.00123

Since the BORGES auxiliary program was not available for authors, a macroscopic cross section approach has been selected for criticality calculations. In this technique, macroscopic cross sections are calculated using the WIMS code with resonance shielding being normally automatically calculated ^[12]. These macroscopic cross sections will be used in the "008" section of CITATION code.

WIMS input files have been prepared for different basic cells including the SFE and CFE fuel plates, Absorber plate, and reflectors (Top, Bottom, and Lateral) based on the six groups of energy structure. Macroscopic cross sections have been extracted from the WIMS output files and implemented in the whole core CITATION model. Important neutronic parameters, including effective multiplication factor as well as forward and adjoint flux distributions, were calculated. Each fuel assembly has been divided into 21 axial mesh points and 5×5 mesh points in the X-Y plane. The first five mesh points and the last five mesh points in axial direction are reserved for top and bottom reflectors, respectively. The reminded 11 mesh points are allocated to active core height.

4.2 Probabilistic approach

In the probabilistic approach, one can rewrite the Eq. (5) as follow ^[8]:

β_{e f f} = \frac{ϕ^{†} {\hat{F}}_{d} Φ}{Φ^{†} \hat{F} ϕ},

(11)

where $\hat{F}$ is the creation operator that takes in to account all neutrons (prompt and delayed) created in the phase space and ${\hat{F}}_{d}$ is the delayed neutron creation operator that takes into account only delayed neutrons. Brackets indicate the integration over the whole phase space. Once again $ϕ$ and $ϕ^{†}$ are forward and adjoint fluxes, respectively, and obtained from eigenvalue equations ^[8]:

\hat{M} ϕ = \frac{1}{k_{e f f}} \hat{F} ϕ,

(12)

\hat{M^{†}} ϕ^{†} = \frac{1}{k_{e f f}} F^{†} ϕ^{†},

(13)

where $\hat{M}$ is loses operator which takes in to account all neutrons leaving the phase space (for example capture or out-scattering) and $\hat{M^{†}}$ has the same meaning for adjoints.

Some techniques have been proposed to evaluate Eq. (11) based on the Monte Carlo method. The most straightforward method is based on the assumption that $≃ ϕ_{p}$ , which means the total flux is approximated by the prompt one. However, this assumption is not crucial since more than 99% of all neutrons are the prompt one. By this assumption, the following eigenvalue equations are obtained instead of Eqs (12) and (13) ^[8]:

\hat{M} ϕ_{p} = \frac{1}{k_{p}} {\hat{F}}_{p} ϕ_{p},

(14)

\hat{M^{†}} ϕ_{p}^{†} = \frac{1}{k_{p}} {\hat{F}}_{p}^{†} ϕ_{p}^{†},

(15)

where ${\hat{F}}_{p}$ is the prompt neutron creation operator. The creation operator $\hat{F}$ can be decomposed into the sum of prompt and delayed ones: $\hat{F} = {\hat{F}}_{p} + {\hat{F}}_{d}$ . By these definitions the Eq. (10) is written as follows ^[8]:

\begin{array}{l} β_{e f f} = \frac{〈 ϕ^{†} {\hat{F}}_{d} Φ 〉}{〈 Φ^{†} \hat{F} ϕ 〉} = \frac{〈 ϕ^{†} (\hat{F} - {\hat{F}}_{p}) Φ 〉}{〈 Φ^{†} \hat{F} ϕ 〉} \\ = 1 - \frac{〈 ϕ^{†} {\hat{F}}_{p} Φ 〉}{〈 Φ^{†} \hat{F} ϕ 〉} ≃ 1 - \frac{〈 ϕ^{†} k_{p} \hat{M} Φ_{p} 〉}{〈 Φ^{†} k_{e f f} \hat{M} Φ_{p} 〉} \\ = 1 - \frac{k_{p} 〈 ϕ^{†} \hat{M} Φ_{p} 〉}{k_{e f f} 〈 Φ^{†} \hat{M} Φ_{p} 〉} = 1 - \frac{k_{p}}{k_{e f f}} . \end{array}

(16)

This technique is known as the "Prompt Method", and its results usually are accepted as a first approximation in Monte Carlo based point kinetic parameters evaluation ^[15]. However, the “Prompt Method” has a weak theoretical justification since it does not use the adjoint weighted parameter in the calculation of $β_{eff}$ ^[8].

5 Results and discussion

At first, cell and core calculations, as well as six groups of energy structure, have been verified via neutronic calculations of the TRR core at different states. The SAR report and MCNP4C probabilistic code have been employed for verification. The deterministic calculation of excess reactivity at each operational state is compared with SAR and cross checked with the MCNP calculation in Table 8. Deterministic excess reactivity calculations in the six groups are in very good agreement with SAR and MCNP results. In the deterministic approach, the maximum relative error is about 6.3% with respect to the SAR report. On the other hand, maximum relative error in the probabilistic approach is only 1.7%. These results confirm input files as well as six groups of energy structure, and they can be used in point kinetic parameters calculation, confidently.

The excess reactivity (pcm) of main states in TRR research reactor, all rods withdrawn

Case	Core State	SAR^[11]	Six GroupWIMS & CITATION(R.E %)	MCNP (R.E %)
1	Cold	6916	6481(6.3)	6934(0.26)
2	0 MW	6549	6398(2.3)	6647(1.5)
3	5 MW (without Xe)	6469	6316(2.4)	6577(1.7)

R.E: Relative Error

Equation (6) has been employed in the deterministic approach to evaluate effective delayed neutron fractions. Forward and adjoint fluxes of each mesh point have been extracted from the CITATION output file and implemented in this equation. As was mentioned earlier, the total delayed neutron fraction ( $β_{eff}$ ) depends on the core isotopic content, especially relative to ²³⁵U and with minor contributions from fission of ²³⁸U. WIMS results showed that 99.5% of all fissions take place in ²³⁵U whilst the fast fission of ²³⁸U is responsible for only 0.5% of fissions.

The MCNP4C code ^[16] is employed in the probabilistic approach to evaluate the effective delayed neutron fraction based on the “Prompt Method” scheme. In this method, the MCNP results of the effective multiplication factor have been calculated in two cases. In the first one, the effective multiplication factor is calculated by considering the TOTNU card of the MCNP code with No as the entry (K_P). Once again, the effective multiplication factor is calculated without the TOTNU card (K_eff). In this regards, the prompt effective delayed neutron fraction is calculated based on Eq. (16).

Calculated values of the effective delayed neutron fractions have been shown in Table 9. The “Prompt Method”, as well as SAR results and experimental data, have been employed for validation.

Effective delayed neutron fractions of TRR reactor core

Groups	Six GroupWIMS & CITATION	Prompt Method	SAR^[11]	Experiment^[17]
1	0.00024	NA	NA	NA
2	0.00158	NA	NA	NA
3	0.00142	NA	NA	NA
4	0.00286	NA	NA	NA
5	0.00085	NA	NA	NA
6	0.00031	NA	NA	NA
Total	0.00726	0.00761	0.00813	0.00786

NA: Not Available

Deterministic, as well as probabilistic, calculations of the total effective delayed neutron fraction are in good agreement with the experiment and SAR report. In comparison with the experimental data, the relative differences of deterministic as well as probabilistic methods are 7.6% and 3.2%, respectively. These quantities are 10.7% and 6.4% respectively in comparison with the SAR report. On the other hand, the deterministic approach has an advantage to provide the delayed neutron fraction in each group beside the total one. The group wise delayed neutron fraction is required in the computational model of some sophisticated dynamic codes, such as PARCS and PARET ^[18].

The mean generation time ( $Λ$ ) and consequently the prompt neutron life time ( $l = Λ k$ ) has been evaluated using Eq. (9). According to this equation, inverse velocity in each energy group is required to calculate the mean generation time. In the WIMS code, the inverse velocity calculation is activated by including the ( $\frac{1}{ν}$ ) absorber (MAT ID 1000) in NREACT and REACT cards. Then, the inverse velocity of each group is calculated according to the following formula:

{(\frac{1}{ϑ})}_{g} = \frac{\sum_{n = 1}^{N m e s h} R R_{a b s o r p t i o n}^{n} V O L_{n}}{2.2 e + 5 \sum_{n = 1}^{N m e s h} ϕ_{g, n} V O L_{n}} .

(17)

The prompt neutron life time of the Tehran Research Reactor has been shown in Table 10. The calculated result is in very good agreement with the SAR one.

Prompt neutron life time of TRR reactor

	Six GroupWIMS & CITATION	SAR
Prompt neutron life time (µ Sec)	48	45

6 Conclusion

Effective point kinetic parameters in the Tehran Research Reactor core have been calculated using deterministic as well as probabilistic schemes. These parameters include the effective delayed neutron fraction and prompt neutron life time. In the deterministic approach, WIMS-D5B and CITATION-LDI2 codes have been implemented efficiently for cell and core calculations, respectively. In order to distinguish between prompt and delayed neutrons, an energy structure in six groups has been selected and group constants have been computed. The six groups' structure was verified via excess reactivity comparison in different states. In the probabilistic approach, the well-known MCNP4C code is employed in the "Prompt Method" scheme. The delayed neutron fraction and mean generation time calculation procedure was verified by SAR, experimental, and “Prompt Method” Monte Carlo results. The results are in very good agreement with each other.

References:

[1]

J. J. Duderstadt, L. J. Hamilton, Nuclear Reactor Analysis, (John Wiley & Sons, Inc, 1976), p. 234

[2]

S. C. Van der Marck, R. K. Meulekamp,

Calculating the effective delayed neutron fraction using Monte Carlo techniques

. Paper presented at PHYSOR 2004- The Physics of Fuel Cycles and Advanced Nuclear Systems: Global Developments Chicago, Illinois, April 25-29, 2004

Baidu Scholar

Google Scholar

[3]

A. Muhammad, M. Iqbal, T. Mahmood, et al.,

Calculation and measurement of kinetic parameters of Pakistan Research Reactor-1 (PARR-1)

. Ann. Nucl. Energy. 38, 44-48 (2011). doi: 10.1016/j.anucene.2010.08.019

Baidu Scholar

Google Scholar

[4]

M. Iqbal, T. Mahmood, Sh. Pervez,

Flow of kinetic parameters in a typical swimming pool type research reactor

. Ann.Nucl.Energy.35, 518-524 (2008). doi: 10.1016/j.anucene.2007.07.024

Baidu Scholar

Google Scholar

[5]

A. Muhammad, M. Iqbal, T. Mahmood,

Kinetic parameters study based on Burn-up for improving the performance of research reactor equilibrium core

. Nucl. Technol. Radiat. 29(4), 253-258 (2014), doi: 10.2298/NTRP1404253M

Baidu Scholar

Google Scholar

[6]

J. E. Hoogenboom,

Methodology of continuous-energy adjoint Monte Carlo for neutron, photon and coupled neutron-photon transport

. Nucl. Sci. Eng. 143, 99-120 (2003), doi: 10.1007/978-3-642-18211-2_98

Baidu Scholar

Google Scholar

[7]

P. K. Meulekamp, S. C. Van der Marck,

Calculating the effective delayed neutron fraction with Monte Carlo

. Nucl. Sci. Eng. 152, 142-148 (2006)

Baidu Scholar

Google Scholar

[8]

V. Becares, S. Perez-Martin, M. Vazquez-Antolin, et al.,

Review and comparison of effective delayed neutron fraction calculation methods with Monte Carlo codes

. Ann. Nucl. Energy. 65, 402-410 (2014), doi: 10.1016/j.anucene.2013.11.038

Baidu Scholar

Google Scholar

[9]

W. M. Stacey, Nuclear Reactor Physics, 2nd Ed. (Wiely-VCH Verlag Gmbh & Co. KGaA, 2007) pp. 145-146

[10]

A. F. Henry,

Nuclear Reactor Analysis

. (MIT University of Technology, 1975)

Baidu Scholar

Google Scholar

[11]

AEOI,

Safety Analysis Report for Tehran Research Reactor

. (2001)

Baidu Scholar

Google Scholar

[12]

NEA, WIMSD-5B, A Neutronic Code for Standard Lattice Physics Analysis. (Code Manual 2003)

[13]

T. B. Fowler, D. R. Vondy, G. W. Cunninghan,

Nuclear reactor core analysis code-CITATION. (USAEC, Report ORNL-TM-2496, ORNL-4078, ReV.2

, 1971)

Baidu Scholar

Google Scholar

[14]

ANL-5800, Reactor Physics Constants. (Argonne National Laboratory, 1963)

[15]

M. Carta, S. Dulla, V. Peluso, et al.,

Calculation of the Effective Delayed Neutron Fraction by Deterministic and Monte Carlo Methods

. SCI. TECHNOL. NUCL. INS. 2011, Article ID 584256, 8 pages (2011)doi: 10.1155/2011/584256

Baidu Scholar

Google Scholar

[16]

J. F. Briesmeister, MCNP4C manual, Monte Carlo N-Particle Transport Code System, RSICC code package CCC-700/MCNP4C, (Oak Ridge National Lab., TN, and DOE, USA 2000.)

[17]

M. Zaker,

Effective delayed neutron fraction and prompt neutron lifetime of Tehran research reactor

. Ann. Nucl. Energy. 30, 1591-1596 (2003), doi: 10.1016/S0306-4549(03)00081-1

Baidu Scholar

Google Scholar

[18]

T. Downar, Y. Xu, T. Kozlowski, et al., PARCS v2.7 U.S.NRC Core Neutronics Simulator User Manual, (Purdue University and RES/U.S. NRC, USA, 2006)