Report on Thermal Neutron Diffusion Length Measurement in Reactor Grade Graphite Using MCNP and COMSOL Multiphysics

0.96 MB·English

by  S. R. Mirfayzi

Report on Thermal Neutron Diffusion Length
Measurement in Reactor Grade Graphite Using
MCNP and COMSOL Multiphysics
3
1
S. R. Mirfayzi
0
2 School of Physics and Astronomy University of Birmingham Birmingham B15 2TT United
Kingdom
n
a E-mail: [email protected]
J

8
Abstract. Neutrondiffusionlengthinreactorgradegraphiteismeasuredbothexperimentally

] and theoretically. The experimental work includes Monte Carlo (MC) coding using ’MCNP’
x and Finite Element Analysis (FEA) coding suing ’COMSOL Multiphysics’ and Matlab.
e
The MCNP code is adopted to simulate the thermal neutron diffusion length in a reactor
-
l moderator of 2m x 2m with slightly enriched uranium (235U), accompanied with a model
c
designed for thermal hydraulic analysis using point kinetic equations, based on partial and
u
ordinarydifferentialequation. Thetheoreticalworkincludesnumericalapproximationmethods
n
[ including transcendental technique to illustrate the iteration process with the FEA method.
Finally collision density of thermal neutron in graphite is measured, also specific heat relation

1 dependability of collision density is also calculated theoretically, the thermal neutron diffusion
v
lengthingraphiteisevaluatedat50.85±0.3cmusingCOMSOLMultiphysicsand50.95±0.5cm
9
usingMCNP.Finallythetotalneutroncross-sectionisderivedusingFEAinaninverseiteration
9
form.
6
1
.
1
0
3
1. Introduction
1
: This work demonstrates an analytic approach accompanied with models of Finite Element
v
Analysis (FEA) and Monte Carlo (MC) with an experimental measure on neutron cross-section
i
X
and slowing down process. In MC approach Monte Carlo N-Particle Transport Code (MCNP)
r is used to simulate the simplified version of reactor moderation process. Similarly in FEA the
a
moderator modelled (Assuming a symmetrical distribution) using point kinetic equations, based
on partial and ordinary differential equation in software package.
2. Theoretical Calculations
Having the number of particles found in a volume element dr where dr = dxdydz at r with a
vector with solid angle dΩ at Ω be donated by [1]:
N(r,Ω,t)drdΩ (1)
Therefore can have:
dN (cid:90)
= −NVσ+ N(r,Ω(cid:48),t)Vσ f(Ω.Ω(cid:48))dΩ(cid:48)+S(r,Ω,t) (2)
s
dt
Where the first term dN donates the number of particles present in given volume (particle den-
dt
sity) and second term −NVσ represents the total number of particles removed from the given
volume by scattering and capture. σ is representing the total cross-section. The third term
represents the total number of particles scattered into the given volume, and f(Ω.Ω(cid:48)) represents
the relative probability of scattering through an angle whose cosine is Ω.Ω(cid:48), where Ω(cid:48) is a unit
vector in the direction of the initial velocity and Ω is unit vector in the final direction. Finally
S(r,Ω,t) is the external source term available in the system and is given by:
N(r,Ω,t)drdΩ = N(Z,ϕ)dZdϕdφ (3)
Where ϕ is the cosine of the velocity vector in the Z direction and φ is the longitude of velocity
vector.
Figure 1. The Velocity Vector: Where ϕ is the cosine of the velocity vector in the Z direction
and φ is the longitude of velocity vector and Ω is unit vector in the final direction.
Now an assumption can be made such:
(cid:90) +1
N (Z) = 2π N(Z,ϕ)dϕ (4)
0
−1
Rewrite the equation 1 as:
δN(Z,ϕ) (cid:90)
V = −N(Z,ϕ)Vσ+ N(Z,ϕ(cid:48))Vσ f(ϕ )dΩ+S(z) (5)
ϕ s 0
dZ
Where ϕ is the cosine of the angle between initial and final velocities and it can be found by:
0
cosθcosθ(cid:48)+sinθsinθ(cid:48)cos(φ−φ(cid:48)) (6)
2
Or using:
(cid:113) (cid:113)
ϕ = ϕϕ(cid:48)+ 1−ϕ2 1−ϕ(cid:48)2cos(φ−φ(cid:48)) (7)
0
Now if the collision function of F(ϕ ) expanded in spherical harmonics:
0
(cid:88)∞ 2l+1
F(ϕ ) = F P (ϕ ) (8)
0 1 1 0
4π
0
(cid:82)
With F = f(ϕ )P (ϕ )dΩ. Using similar expansion the phase density function will be:
1 0 1 0
(cid:88)∞ 2l+1
N(Z,ϕ) = N (Z)P (ϕ) (9)
1 1
4π
0
(cid:82)
Where N(Z,ϕ) = N(Z,ϕ)P (ϕ)dΩ, with assumption that N(Z,ϕ) is isotropic, three condi-
1
tions must be satisfied all the times: first, it is far from the source (equal to Mean Free Path
(MFP)),second,itisfarfromtheboundaries; third,theprobabilityofcaptureissmallcompared
to probability of scattering. Having all the conditions satisfied, the following can be assumed:
1
∼
N(Z,ϕ) = (N (Z)+3ϕN (Z)) (10)
0 1
4π
(cid:82)
Where the second term in the bracket donates the particle flux (J). Here N = ϕN(Z,ϕ)dΩ =
1
J/V. For simplicity we choose our unit such that V = 1 and σ = 1, hence:
σ
s
1−f = (11)
σ
Now the Boltzmann equation takes the form of:
dN (cid:90)
N = −N +(1−f) N(Z,ϕ(cid:48))f(N )dΩ(cid:48)+S(Z) (12)
0
dZ
By integrating the equation over all possible angles (dΩ) we have:
dN
1
= −N +(1−f)N +4πS(Z) (13)
0 0
dZ
f is normalized in such a way that (cid:82) f(Ω.Ω(cid:48))dΩ(cid:48) = (cid:82) f(Ω.Ω)dΩ = 1, hence going back to Eq. 9
for the case l = 0 we have:
(cid:90)
F = f(ϕ )dΩ = 1 (14)
0 0
Hence by integrating over all angles and Multiplying by ϕ we have:
1dN
0
= −N +(1−f)F N (15)
1 1 1
3 dZ
3
Now the second order differential equation gives:
1 d2ϕ
0
− = −fN +S (Z) (16)
3(1−f dZ2 0 0
1
This also can be written as:
1 1
(cid:53)2N − N + S = 0 (17)
0 L2 0 D 0
Eqs. 16 and 17 are known as diffusion equation. Here L is diffusion Length abd D is diffusion
coefficient and it is equal to 1 λs = λtrV, where λ and λ are the scattering and transport
31−f1 3 s tr
mean free path. λ can be calculated from:
tr
λ
s
λ = (18)
tr
1−(cosθ)
av
and (cosθ) is equal to (cid:82) f(ϕ )ϕ dΩ = f . Also L2 can be measured using following relation:
av 0 0 1
λ λ
L2 = c tr (19)
3
Here λ is the capture mean free path.
c
MAXIMUM ENERGY LOSS
If a neutron with initial velocity V collides with a nucleus of mass M (at rest), then in the
0
Centre of Mass (CoM) system, the initial velocity is MV /M+1 after collision. The momentum
0
of of neutron and the nucleus will be equal to magnitude oppositely directed vector. Figure 2
demonstrates the collision in CoM system.
As demonstrated in Fig. 2 the θ is the deflection angle and Θ is angle on the final velocity v.
The v2 in this case is given by:
Mv v
0 0
cosθ+ = vcosΘ (20)
M +1 M +1
Mv v 2Mv2
( 0 )2+( 0 )2− 0 cosθ = v2 (21)
M +1 M +1 M +1
so
(M +1)2 v
cosθ = 1− 1− )2 (22)
2M v
0
since u = logE0 then:
E
(M +1)2
cosθ = 1− 1−e−u (23)
2M
4
Figure 2. The Collision in Centre-of-Mass System:If a neutron with initial velocity V collides
0
with a nucleus of mass M at rest, then in the Centre of Mass (CoM)system the initial velocity
is MV /M +1 after collision. The momentum of of neutron and the nucleus will be equal to
0
magnitude oppositely directed vector.Here θ is the deflection angle and Θ is angle on the final
velocity v.
now differential cross-section gives:
dcosθ (M +1)2
= − e−u (24)
du 2M
Hence:
(M +1)2 M −1
−u u
cosΘ = − e 2 − e2 (25)
2 2
Therefore the maximum logarithmic energy loss can be calculated from:
M +1
q = log( )2 (26)
M
M −1
The q is at most when Θ = π. Now going back to the problem we can redefine the collision
M
density function as:
(M +1)2 (M +1) M −1
F(ϕ0,u) = e−u×δ(ϕ0−( e−2u − eu2) (27)
8πM 2 2
The term (M+1)2e−u is the normalization constant chosen to satisfy (cid:82) dΩ(cid:82) duf(ϕ ,u) = 1 and
8πM 0
5
(cid:82)
δ is the Dirac δ function. So that δ(x−a) = 0 when x (cid:54)= 0 and d(x−a)F(x)dx = F(a). Now
the average logarithmic loss (ξ) can be calculated from:
(M +1)2
ξ = 1− q e−qm (28)
m
4M
and (cos) v is 2/3M.
a
Energy Distribution of Slowed Down Neutrons
I. Stationary Case
The average collision density per unit time, with logarithmic energy intervals is given by:
(cid:90) u
ψ (u) = du(cid:48)ψ (u(cid:48))h(u(cid:48))f (u−u(cid:48))+δ(u) (29)
0 0 0
0
where f (u) is (M +1)2e−u/4m for u ≤ q and it is zero otherwise. In stationary case the total
0 m
number of neutron produced is unity per unit in this case, i.e. for M = 12, f +0(u) becomes
3.5e−u. Hence the equation 29 becomes:
(cid:90) u
ψ (u) = du(cid:48)ψ (u(cid:48))h(u(cid:48))3.5e−(u−u(cid:48))+δ(u) (30)
0 0
u−qm
where q = 0.72
m
II. Time-dependent Case
The time dependent when there is no absorption in the system and source strength is unity and
is given by:
l(u)dψ (cid:90) u
0 +ψ (u,t) = du(cid:48)ψ (u(cid:48),t)e−(u−u(cid:48))+δ(u)δ(t) (31)
0 0
v δt
0
where l(u) is the mean free path and if the mean free path is constant, the Laplacian form of
the equation for M (cid:54)= 1 becomes:
sl (cid:90) u
1+ 0φ (u,s) = du(cid:48)φ (u(cid:48),s)f (u−u(cid:48))+δ(u) (32)
0 0 0
v
0
now:
2 (cid:90) w w(cid:48)φ(w(cid:48),s)
φ(w,s) = dw(cid:48) (33)
(1−r2)w2 (1+w(cid:48))
rw
Eq. 33 applies for u > q where w = l s/v, r = M −1/M +1 and φ(w,s) = (1+w)φ (u,s), so
m 0 0
that the mean free path is proportional to velocity.
III. Rigorous Numerical Solution
The slowing down process is not an easy approach, therefore a more discrete form of solution
also could be defined using:
(cid:90) ∞
F(E) = (cid:88)(E(cid:48) → E)φ(E(cid:48))dE(cid:48)+δ(u) (34)
E
S
6
Where (cid:80) (E(cid:48) → E) is the scattering term between energies E(cid:48) and E. Recalling (cid:80) :
S S
(cid:88)(E(cid:48) → E) = (cid:88)(E(cid:48))P(E(cid:48) → E) (35)
S S
WhereP(E(cid:48) → E)istheprobabilityofcollisionhappensbetweenE(cid:48) andE. Itcanbedefinedby:
P(E(cid:48) → E).E(cid:48)(1−α) = 1 (36)
Hence:
1
P(E(cid:48) → E) = (37)
(1−α)E(cid:48)
Now:
(cid:88)(E(cid:48) → E) = (cid:80)S(E(cid:48)) (38)
(1−α)E(cid:48)
S
For M = 1 the collision is defined as:
(cid:90) ∞ (cid:80) (E(cid:48))
F(E) = S φ(E)dE(cid:48)+δ(E) (39)
E(cid:48)
E
Also the solution with capture process:
F (E) = (cid:80)S(E0)S0exp(−(cid:90) E0 (cid:80)S(E(cid:48))dE(cid:48)) (40)
c (cid:80) (E) E (cid:80) (E(cid:48)) E(cid:48)
t E t
Also For M (cid:54)= 1 the collision density can be found:
(cid:90) E (cid:80) (E(cid:48)) δ(E)
F(E) = S φ(E(cid:48))dE(cid:48)+ (41)
(1−α)E(cid:48) (1−α)E
E/α 0
This only applicable if αE < E < E . A theoretical calculation is performed for an arbitrary
0 0
system and Fig 3 is derived. For graphite the collision density is also measured for different
neutron energy range as demonstrated by Fig. 4.
and for graphite:
The oscillations are due to Plaezack Oscillations which is the fundamental phenomenon
associated with the neutron slowing-down [2]. And finally for the M (cid:54)= 1 with capture:
(cid:88) (cid:88)
F(E) = ( (E)+ (E))φ(E)
s a
(cid:32)(cid:90) E (cid:80) (E(cid:48))φ(E(cid:48)) S (cid:33)
= S dE(cid:48)+ 0 (42)
(1−α)E(cid:48) (1−α)E
E/α 0
(cid:90) E (cid:80) (E(cid:48)) F(E(cid:48)) S
= S dE(cid:48)+ 0
(cid:80) (E(cid:48)) (1−α)E(cid:48) (1−α)E
E/α t 0
7
Figure 3. The Maximum Logarithmic Energy Loss vs. Collision Density
Figure 4. The Collision Density vs. Neutron Energy in Graphite
3. Model Set-up in MCNP
The MCNP code is developed in Los Alamos National Laboratory and it is well-known for
analysing the transport of neutron and γ-rays in matter. MCNP is a continuous energy mod-
ellerwithgeneralizedgeometrytimedependentcodethatimplementsdatafromnuclearlibraries
such as, Evaluated Nuclear Data File (ENDF), Evaluated Nuclear Data Library (ENDL), Acti-
vation Library (ACTL).
The code structure is divided into four main sections. geometry definitions, surface definitions,
material cards, and tallies. Geometry of MCNP is a three-dimensional form defined using cell
and surface cards. For instance Fig 5 demonstrates the geometry setup in this system. Figure 5-
illustrates the reactor moderator, where each cylinder represents the fuel rod containing slightly
enriched 235U. The moderator is reactor grade graphite.
The user can instruct the code to make various analysis with tally cards. The tallies are to
measure the particle current on the surface, particle flux and energy deposition. In fact any
8
Figure 5. The MCNP Geometry Set up: Figure demonstrates the reactor moderator module.
The dimension was set 2m x 2m
quantity in form of Eq. 43 can be tallied [3].
(cid:90)
C = φ(E)f(E)dE (43)
Here φ(E) represent the particle flux and f(E) is the cross-section quantities given in the li-
braries. Table 1 demonstrates the six MCNP standard tallies.
In MCNP when neutron collides with a nucleus: the nuclide will be identified depends on the
Table 1. MCNP Neutron Tallies
Property Data
F1:N Surface Current
F2:N Surface Flux
F4:N Track Length Estimate of Cell Flux
F5a:N Flux at a Point
F6:N Track Length Estimate of Energy Dependence
F7:N Track Length Estimate of Fission Energy Dependence
preferences of target, that is either the S(α,β) treatment or velocity of target; therefore the
nucleus will be sampled for low energy neutrons; neutron capture or absorption will be modelled
9
and either elastic or inelastic reaction depend on the model performance.
However sometimes different nuclide form a material, (where the collision occurs) therefore we
can have:
k−1 n k
(cid:88)(cid:88) (cid:88)(cid:88) (cid:88)(cid:88)
< ξ ≤ (44)
i=1 ti i=1 ti i=1 ti
Where (cid:80) is the microscopic total cross-section of nuclide i. The total cross-section is sum of
ti
the capture cross-sections in the cross-section reference table.
The collision between thermal neutrons and the target will be effected by thermal motion of the
atoms, chemical binding and lattice structure of the target. This is called Free Gas Thermal
Treatment. Hence the effective scattering cross-section in laboratory system is given by [4]:
1 (cid:90) (cid:90) dϕ
σeff(E) = σ (v )v P(v)dv t (45)
s v s ref rel 2
n
Here v is particle velocity, v is the relative velocity, P(v) is the probability density function
n rel
and ϕ as explained before is the cosine angle of velocity vector. The relative velocity can there-
fore is given by:
v = (v2 +v2−2v vϕ )1/2 (46)
rel n n t
The density function is also given by:
4
P(v) = β3v2e−β2v2 (47)
π1/2
where β = (AMn)1/2. However most of the time in equation 45 the σ can be ignored for heavy
2kT s
nuclei, where σ can have moderating effect and is given by [4]:
rel
(cid:113)
P(v,ϕ) ∝ v2v2−2v.v ϕ v2e−β2v2 (48)
n n t
In MCNP there are also two types of capture, analogue and implicit. Analogue occurs when
the particle is killed with probability of σa. Where σ and σ is the absorption and total cross-
σt a t
section respectively. Implicit capture happens when neutron weight (W ) is reduced by number
n
of collisions and is given by:
σ
a
W = (1− )W (49)
n+1 n
σ
t
The elastic scattering directed by two body kinematics:
1
E = E n((1−α)Θ m+1+α) (50)
out i c
2
Where Θ m is the center of mass cosine of angle between incident and existing path direction.
c
Where in inelastic an scatter the particle reaction is chosen such as (n,n(cid:48)), (n,2n), (n,f), and
(n,n(cid:48)α), and is given by [4]:
E +2Θ(A+1)(cid:112)EE(cid:48)
E(cid:48) = E(cid:48)m+( cm) (51)
c (A+1)2
10