Prompt Time Constants of a Reflected Reactor*Tao Ye

Chaobin Chen Weili Sun Benai Zhang

Dongfeng Tian

Institute of Applied Physics and Computational Mathematics, Beijing, 100088, China*E-mail: yetao@aees.kyushu-u.ac.jp

Based on G. D. Spriggs’ two-region kinetics model, a two-group point reactor kinetics model isdeveloped. With the help of MCNP code, the modified model calculates prompt time decay constants of one
benchmark reactor, PU-MET-FAST-024. The results of fundamental and secondary modes agree well withMCNP time fitting results in different subcritical reactivities.

1. INTRODUCTION
Time eigenvalue of transportation equation, alpha, is defined to describe all neutrons’ time behavior(increasing or decreasing) in a nuclear reactor. Its number reflects the criticality also. The time constant,especially prompt time constant, had been studied for 60 years. Lots of reflected reactor’s experimental
data cannot be satisfactorily explained using the standard point kinetic model 1. And multiple decay modes
near delayed critical were also observed, which of course cannot be described by standard point kineticmodel.
The existing numerical transportation codes, such as MCNP4C 2 and TART 3, 4 , can do the job wellwith only the fundamental mode calculated. By using alpha static criticality method, MCNP4C is a goodtool if keff is close 1, which means the reactor is near delayed critical. But MCNP4C’s calculation may bevery difficult and time-consuming if the reactor has more negative reactivity or reflector contains hydrogen,or both.
In the region of analytical method, many works contains too much mathematics, which are not easy tocalculate and compare with experimental data. G. D. Spriggs’ one-group, two-region kinetic model basedon Avery-Cohn model is simple, calculable. The model introduces simple probability relationships essential
to calculating the coupling parameters between core and reflector,1 and derives the reflected-core inhour
equation which contains multiple decay modes. However, Spriggs model cannot well describe multiple
time constants of the thermal reflected reactor. In this kind of reactor, thermal neutrons with long lifetime
contribute much to the time constant. Because of importance of thermal neutrons in such fast-thermalreactor, we present a simplified two-group, two-region kinetic model (2G2R) based on Spriggs model, andrewrite the reflected-core inhour equation. With the help of MCNP code, we calculated the coupling
parameters, neutron lifetimes and first and secondary time constant of a spherical benchmark reactor,
PU-MET-FAST-024.6 Because we don’t have experimental data, the results of time constants are also
compare with 3 different models, MCNP time fitting method, alpha static method (MCNP4C), and Spriggsmodel. The results of 2G2R model agree well with MCNP time fitting method which can be thought as anexperiment in computer.

2006 Symposium on Nuclear Data

1
SND2006-V.09-1

2. ALPHA STATIC CRITICALITY METHODMCNP4C code introduced a new feature to calculate the fundamental mode of prompt time
eigenvalue2. It is based on alpha static criticality method. In subcritical condition, the equation is⎛υΩ ⋅ ∇N + ⎜⎜ σ t +⎝⎞⎛α⎞α⎟ ⋅ υN = ∫∫ dΩ' dE ' ⎜ σ s 'υ ' N '+2 ⋅ υ ' N ' δ (Ω'−Ω)δ ( E '− E ) ⎟⎟⎜⎟υ ⎠υ⎠⎝.

(1)⎛⎞α⎞α⎟ ⋅ υN = ∫∫ dΩ' dE ' ⎜ σ s 'υ ' N '+2 ⋅ υ ' N ' δ (Ω'−Ω)δ ( E '− E ) ⎟⎜⎟υ ⎟⎠υ
⎝
⎠.

(1’)

1+ ∫∫ dΩ' dE ' (vσ f υ ' N ')k'The calculation procedures are to get k’≈1 by searching proper alpha. Then equation becomes⎛υΩ ⋅ ∇N + ⎜⎜ σ t +⎝+ ∫∫ dΩ' dE ' (vσ f υ ' N ')which is the alpha eigenequation.
If a reactor has more negative reactivity or reflector contains hydrogen, or both, the ratio between (α/v)
term and σt term can be very large, which will results non-physical high particle weight and stops thecalculation. We add an adjusting parameter to lower the ratio. The modified equation is⎛υΩ ⋅ ∇N + ⎜⎜ σ t + η⎝⎞⎛α⎞α⎟ ⋅ υN = ∫∫ dΩ' dE ' ⎜ σ s 'υ ' N '+(1 + η ) ⋅ υ ' N ' δ (Ω'−Ω)δ ( E '− E ) ⎟⎟⎜υ ⎟⎠υ⎠⎝1+ ∫∫ dΩ' dE ' (vσ f υ ' N ')k'(2)The modification can only weaken the tendency of getting huge particle weight, and restrictedly
extend the usage of MCNP4C. The k’ intends to converge to a number larger than unit if reactor is in a
deeper subcriticality, which means calculated alpha is smaller than true value in number axis. The
determination of adjusting parameter is a little arbitrary based on various calculation conditions. Once
confirmed, it shall not change in the running.

3. MCNP TIME FITTING
In time dependant transportation equation, neutron density or flux has the formal solution,5()()

r rr rα t
N r , Ω, E , t = ∑ N j r , Ω, E × e j .
j =0(3)

In a subcritical system, all αj values are negative. We assume αj‘s absolute values increase with increasing j.
And α0 is the largest one, the fundamental time constant. If we add a pulse source at zero time, theneutron’s time distribution will start a buildup in the beginning, then drop to multiple decay mode which
has nothing to do with source anymore.
By integrating volume, solid angle and energy, the current term becomes leakage term. But its time
behavior still follows formula (3). which means we can use MCNP’s tally option, F1, to count leakage
neutrons’ time distribution as system’s time distribution. Then, we use formula (3) to fit time distribution toget multiple time constants. The fitting coefficients are time-independent flux corresponding to each decay
mode. However, the coefficient, Nj, is not concerned in this work.
With enough neutron source particles (NPS) and adequate time, MCNP F1 tally can explain Rossi-α
measurement well in the vicinity of delayed critical. At this point, MCNP F1 tally is doing the same thing

2006 Symposium on Nuclear Data

2SND2006-V.09-2

as Rossi-α does. The difference is Rossi-α has background term, which is not easy to get rid of fromexperimental counts, and will conceal the lower decay modes, especially α0 , in a deeper subcritical systemwith a thermal reflector. To the contrary, MCNP F1 tally can display all lower decay modes without
interference of background term. So, MCNP F1 tally can be seen as an imaginary Rossi-α measurement incomputer.

4. 2G2R MODEL
We adopt the conventional diffusion approximation to deal with a two-region system consisting of acore surrounded by a non-multiplying, source-free reflector.1 The simplified model can be described as thefollowing set of two-group coupled differential equations.

⎛ dφ1c ⎞⎟⎜⎜ υ1c dt ⎟⎜ dφ ⎟ ⎛ − Bc2 D1c + (1 − β eff )vΣ f − Σ r
1c1c⎜ 2c ⎟ ⎜Σ s12 c⎜ υ 2 c dt ⎟ ⎜2
⎜ dφ ⎟ = ⎜Bc D1c f cr11
⎜ 1r ⎟ ⎜0⎜ υ1r dt ⎟ ⎜⎝⎜ dφ ⎟⎜⎜ 2 r ⎟⎟⎝ υ 2 r dt ⎠(1 − β )vΣeff

Br2 D1r f rc11f2c− B D2c − Σ a2 c02
Bc D2c f cr222c⎞⎛ φ1c ⎞⎟⎜ ⎟B D2 r f rc22 ⎟⎜ φ2c ⎟⎟⎜ ⎟0⎟⎜ φ1r ⎟2− Br D2 r − Σ a2 r ⎟⎠⎜⎝ φ2 r ⎟⎠02r

0− B D1r − Σ r1r
Σ s12 r
2r

(4)
Where subscript 1 represents the fast group (E>1eV), and 2 represents thermal group (E<1eV). c means
core, and r means reflector. For simplification, we only include effective fraction of delayed neutrons, βeff ,in set of equations.
After Laplace transformation, we can get inhour equation,[(ω l + 1 − k (1 − β ))(ω l + 1) − f ]× [(ω l + 1)(ω l + 1) − f ].
− k (1 − β )[(ω l + 1)(ω l + 1) k + f ] = 0
1c1c

eff

1r

2c

112c2r

eff

1r

2r

22s12 c(5)12In most reflected systems, the thermal neutron’s number is a few orders smaller than fast neutron’s number,which leads to a neglectable f22(≈0). And average thermal neutron’s lifetime is sufficiently small such thatωj l2c<<1 for all possible j roots. Introducing the definition of reactivity, the inhour equation can berewritten asω

f12 k 2 c (1 − β eff ) ⎡⎤l1c
1 − f11
f11 l1r1−ω
−= ρ − β eff⎢1 −⎥ = 1 − β eff −(ω l1r + 1)(kc + Δ )kc + Δkc + Δkc + Δ⎣ (ω l1r + 1)(ω l 2 r + 1) ⎦(6)
In many cases, f22 and ks12c can be neglected and l2r is very large such that ωj l2r>>1, then the inhour
equation is back to Spriggs model shape with one feedback constant, and only describe fast neutron’s timeconstant. But we will lose a root related with l2r by this simplification. It is better to resolve cubic equation[(ω l + 1 − k (1 − β ))(ω l + 1) − f ]× (ω l + 1) − k (1 − β ) f = 0 .
1c1ceff

1r

112r2ceff12(7)with determined coefficients calculated by MCNP code.

5. MODEL CALCULATIONS

2006 Symposium on Nuclear Data

3
SND2006-V.09-3

Because we don’t have experimental data to test the 2G2R model, a benchmark model,PU-MET-FAST-024,6 is chosen as preliminary test. A summary of the reactor is given in Table I. And 3
different models do the same calculations for comparison between them. One model, MCNP time fitting, isregarded as imaginary experiment because of similarity with Rossi-α measurement. All results are listedbelow.
Table I: Simplified PUT-MET-FAST-024 Reactor Description 6Region/Dimension

Material

Core (spherical)
6cm radius

Reflector(spherical shell)
1.55cm thickAtom Density(×10-24) cm-3239

Pu

3.6620×10-2240

Pu6.6944×10-4

Ga

2.1962×10-3

Fe

1.4126×104O

2.8972×10-4

Ni1.9748×10-3

C

3.8814×10-2

H

7.7616×10-2

D

1.1644×10-5

Table II: Integral Quantities with different reflector thick
Model / Reflector’s thick

0cm
Bare Reactor

0.6cm

1.0cm

1.55cm

0.95323

0.97277

0.99823keff
MCNP

0.92311

Fundamental Time Constant, α0 (μs-1)
Alpha Static Method-25.88-2.372-0.45-0.035MCNP Time Fitting-29.48-1.23-0.53-0.038

Spriggs Modela-22.43-10.0-3.68-0.16

2G2R Modela-22.40

-0.57

-0.26

-0.20

-1

Secondary Time Constant, α1 (μs )MCNP Time Fitting－-8.26

-5.10

-0.37

2G2R Model－-10.21

-4.55

-0.82

a－βeff=0.0023 is taken.

Figures a to d display the MCNP time fitting results. Fig. a is for bare system, and Fig. d is for criticalsystem. In Fig. b and c, neutron counts shows a sharp drop of 2 orders in a few microseconds and still do
not enter the fundamental decay mode. Rossi-α measurement will be difficult to get the fundamental modefor such reactors, because secondary decay mode is dominant and background neutrons cover and destroy
the fundamental decay mode.

2006 Symposium on Nuclear Data

4SND2006-V.09-4

101

1

Neutron Counts (μs )

-1

0.1

Reflector Thick is 0.6cm.
MCNP F1 Tally
Two Decay Mode Fitting
0.4 ~ 10 (μs) region-1
α0= -1.23 (μs )

0.1

-1

Neutron Counts (μs )100.01

1E-3

Reflector Thick is 0cm.
MCNP F1 Tally
One Decay Mode Fitting
0.11 ~ 0.5 (μs) region-1
α0= -29.48 (μs )1E-4

1E-50.01-1α1= -8.26 (μs )

1E-3
1E-4
1E-5
1E-6
1E-7

1E-60.0

0.1

0.2

0.3t (μs)

0.4

00.5

2468t (μs)

a1012b10

1

Reflector Thick is 1.55cm.
MCNP F1 Tally
Two Decay Mode Fitting
1.1 ~ 350 (μs) region-1
α0= -0.038 (μs )0.1

Neutron Counts (μs )-1

Neutron Counts (μs )

1-1

0.01
1E-3

Reflector Thick is 1.0cm.
MCNP F1 Tally
Three Decay Mode Fitting
0.01 ~ 12 (μs) region-1
α0= -0.53 (μs )

1E-4
1E-5
1E-60.1-1

α1= -0.37 (μs )0.01-1

α1= -5.10 (μs )

1E-7
1E-8

1E-30246810t (μs)

c12141618200102030405060708090100t (μs)

dFig. a. Leakage neutron time distribution and its fitting with 0cm-thick reflector.
Fig. b. Leakage neutron time distribution and its fitting with 0.6cm-thick reflector.
Fig. c. Leakage neutron time distribution and its fitting with 1.0cm-thick reflector.
Fig. d. Leakage neutron time distribution and its fitting with 1.55cm-thick reflector.

6. DISCUSSIONComparing the results given in Table II, we can see that all 2 time constants of 2G2R model are close
to MCNP time fitting method, the imaginary experiments. Two reasons contribute much to this. Firstly,
two-group calculation is included in 2G2R model. Second, all coupling parameters from Spriggs model are
determined by MCNP running.
At the same time, fundamental time constants calculated by Spriggs model results are close to MCNP
time fitting for critical and bare system, and are not for two systems in the middle. Two middle systems’
results are close to MCNP time fitting’s secondary time constant. These two features can be explained thatSpriggs model’s time constant reflect the dominant time decay behavior which may not be the fundamental
decay mode.
Though alpha static method (MCNP4C) ‘s results are close to MCNP time fitting, the convergence ofk’ in alpha static equation (1) or (2) is departure from unit with a few percent error for two middle systems,
which strongly lowers the results of alpha static method.

2006 Symposium on Nuclear Data5
SND2006-V.09-5

According to the discussion above, 2G2R model provides a simple way to analyze multiple time decay
modes quantitatively.

REFERENCES
[1].
[2].
[3].
[4].
[5].
[6].

G. D. Spriggs et al., “Two-Region Kinetic Model for Reflected Reactor”, Ann. Nucl. Energy, 24,
205, (1997).
Judith F. Briesmeister, “MCNPTM-A General Monte Carlo N – Particle Transport Code”,
LA-13709-M, Arpi. 10, (2001).
Dermott E. Cullen, “TART 2002: A Coupled Neutron-Photon 3-D Combinatorial Geometry, TimeDependent Monte Carlo Transport Code”, UCRL-ID-126455, Rev. 4, June (2003).
Dermott E. Cullen et al., “Static and Dynamic Criticality: Are They Different?”,
UCRL-TR-201506, Nov. 22, (2003).
G. I. Bell, S. Glasstone, Nuclear Reactor Theory, Van Nostrand, Princeton, N. J., (1970).
M. V. Gorbatenko et al., “Polyethylene-Reflected Spherical Assembly of 239PU(δ, 98%)”, Vol. I,
NEA/NSC/DOC/(95).

APPENDEX Effective multiplication factor, keff
To calculate effective multiplication factor, keff , we resolve the set of equations below,− Bc2 D1c +1(1 − β eff )vΣ f1c − Σ r1ck eff
Σ s12 c1(1 − β eff )vΣ f2 ck eff− Bc2 D2 c − Σ a2 c

Br2 D1r f rc1100

Br2 D2 r f rc22

0− Br2 D1r − Σ r1r

0B D2 c f cr22

Σ s12 r− B D2 r − Σ a2 r

Bc2 D1c f cr112c0= 0 (A1)2r

The effective multiplication factor, keff and kc are,k eff =vΣ f 2 c Σ a2 c
⎛ Σ s12 c Σ r1c⎞⎤1 ⎡ vΣ f1c Σ r1c⎜⎟⎥++f⎢12⎟(1 − f11 ) ⎣⎢ 1 + τ 1 Bc2 1 + τ 2 Bc2 (1 − f 22 ) ⎜⎝ 1 + τ 1 Bc2⎠⎦⎥(())

k 2 c k s12 c + f12 ⎤ k c + Δ⎡1=≈
⎢k1c +(1 − f11 ) ⎣(1 − f 22 ) ⎥⎦ (1 − f11 )k c = k1c + k 2 c k s12 c =vΣ f1c Σ r1c

1 + τ 1B

2006 Symposium on Nuclear Data

2c+vΣ f 2 c Σ a2 c Σ s12 c Σ r1c

1 + τ 2 Bc2

1 + τ 1 Bc2

.

(A2)

(A3)

6
SND2006-V.09-6