Download Neutron Diffusion: Thermal and Fast Neutron Values and Maxwell Distribution and more Study notes Engineering in PDF only on Docsity!
PHGN 590
Introduction to Nuclear Reactor Physics
J.A. McNeil
Modeling Neutron Diffusion in Reactors
(February 29, 2009)
Packages to be Used
Data
Cons = 8 k B Ø 1.38066 μ 10 ^ - 23 , T room Ø 293.15, e -> 1.60219 μ 10 ^ - 19 , m n Ø 1.674929 μ 10 ^ - 27 < ;
D2Odata = 8 r Ø .001105, n d -> .03323, Ss -> .4519, Sg Ø 4.42 μ 10 ^ - 5 , Sf Ø 0 , n Ø 0 < ;
C12data = 8 r Ø .00160, n d Ø .08023, Ss -> .3811, Sg Ø .0002728 , Sf Ø 0 , n Ø 0 < ;
H* Thermal neutron values *L
Nadata = 8 r Ø .00097, n d Ø .02541, Ss Ø .08131, Sg Ø .01347 , Sf Ø 0 , n Ø 0 < ;
U235data = 8 r Ø .01886, n d Ø .04833, Ss Ø .01588, Sg Ø 4.833, Sf -> 28.37, n Ø 2.42 < ;
U238data = 8 r Ø .0191, n d Ø .04833, Ss Ø .4301, Sg Ø .13194, Sf -> 0 , n Ø 0 < ;
Pu239data = 8 r Ø .0196, n d Ø .04938, Ss Ø .3902, Sg Ø 13.27, Sf -> 36.66 , n Ø 2.98 < ;
H* Fast neutron values *L
Nadata = 8 r Ø .00097, n d Ø .02541, Ss Ø .083853, Sg Ø .000020328 , Sf Ø 0 , n Ø 0 < ;
U235data = 8 r Ø .01886, n d Ø .04833, Ss Ø .328644, Sg Ø .0120825, Sf -> .06766, n Ø 2.6 < ;
U238data = 8 r Ø .0191, n d Ø .04833, Ss Ø .33347, Sg Ø .007732, Sf -> .004591, n Ø 2.6 < ;
Pu239data = 8 r Ø .0196, n d Ø .04938, Ss Ø .33578, Sg Ø .0128388, Sf -> .091353, n Ø 2.98 < ;
Introduction
Nuclear energy arises from the splitting of large nuclei by neutrons which, in addition to releasing copius amounts of energy,
release more neutrons. These in turn can cause additional fissions leading to a sustained chain reaction, or in the case of a nuclear
weapon, an exponentially increasing cascade and explosion.
Fundamentally, a neutron interacting with a nucleus can undergo three fundamental scattering processes: 1) it can bounce off
(elastically or inelastically), 2) it can be absorbed by emitting a gamma, and 3) for fissile nuclei, it can be absorbed and then
induce fission which releases additional neutrons. In essence the control of nuclear energy comes down to controlling the
neutrons. The relative probabilities for any of these events are governed by the cross section. The rate for any process is the
neutron flux (number density times average speed) times the relevant cross section times the number density of target nuclei. For
convenience, these last two factors are usually multiplied together to give the "macroscopic cross section". Specifically, the rate
that neutrons are absorbed at location r
is given by Sa j Hr
L,where j Hr
Lis the neutron flux (number of neutrons per unit area
per unit time) and Sa is the macroscopic cross section for absorption. (For fissile materials the absorption cross section includes
both the gamma and fission processes.) Using Fick's law and the continuity equation, one arives at the governing equation for the
neutron flux,
D “
”^2
j Ir
, tM + Hn Sf - SaL j Ir
, tM =
v
ê
∂ j Ir
, tM
∂ t
where j is the neutron flux (neutron number density times average speed), D is the diffusion constant due to elastic scattering, n is
the average number of neutrons produced by a fission, Sa and Sf are the macroscopic cross sections for absorption and fission
respectively. The diffusion constant is related to the macroscopic scattering cross section by
D = 1 ê H 3 HSa + Ss H 1 - mLLL.
(This is an approximate expression.) Analysis of the problem can be simplified by introducing a parameter, k, which scales the
rate of neutron production. One can consider k to be adjusted to give a steady state time-independent neutron flux,
D “
j Hr
L +
n Sf
k
- Sa j Hr
L = 0.
As will be explained later, k is the net neutron multiplication factor that determines if the reactor is critical (self-sustaining).
Define the "buckling" constant,
B
D
n Sf
k
- Sa ,
and, solving for k find,
k =
n Sf
D B
+ Sa
In terms of the buckling constant, the steady-state diffusion equation can be written as an eigenvalue problem,
”^2
j Hr
L = - B
j Hr
L.
One strategy to investigate criticality of a reactor configuration is to solve the eigenvalue problem to obtain allowed values of the
buckling constant, B, for some reactor composition and geometry. The geometric information about the reactor is contained in B.
From Eq.(5) one can understand that the lowest buckling eigenvalue will dominate any cascade process since any larger eigen-
value necessarily yields a smaller neutron multiplication factor. Setting k = 1 for the lowest eigenvalue will determine the
geometry that will give a critical reactor for the given fuel/moderator mixture.
Maxwell velocity distribution
Maxwell distribution normalization check = 1
Vavg =
2
p
mn
T kB
Vrms = 3
T kB
mn
Vpeak =
2 T kB
mn
v
P
Maxwell Veclocity Distribution
Clear @ A, v, vboost, vcm, En D
vpeakvalue = 100 vpeak @ T roomD ê. Cons; H* most probable speed in cm ê s *L
vrmsvalue = 100 vrms @ T roomD ê. Cons ; H* rms speed in cm ê s *L
vavgvalue = 100 vavg @ T roomD ê. Cons; H* average speed in cm ê s *L
vspeed = vavgvalue;
vboost @ A _D = vspeed ê H 1 + A L 8 0 , 0 , 1 < ; H* speed between CM and lab frames *L
vcm @ A _D = A vspeed ê H 1 + A L ;
En = H k B T room ê e L ê. Cons;
Print @ " Thermal energy, k B T = ", En, " eV" D
Print @ " V avg = ", vavgvalue, " cm ê s" D ; Print @ " V rms = ", vrmsvalue, " cm ê s" D ;
Print @ " V peak = ", vpeakvalue, " cm ê s" D ;
H* v = 2.2 10 ^ 5 cm ê s *L
Thermal energy, kB T = 0.0252617 eV
Vavg = 248 062. cmês
Vrms = 269 247. cmês
Vpeak = 219 839. cmês
Task 1: Point Source
Task 1. Before tackling nuclear reactor models, let's just study simple neutron diffusion without any fission. Consider an infinite
volume of some material, say sodium to be specific, at room temperature with a thermal neutron point source at the origin that
emits S neutrons per second isotropically. The analytic solution for the steady state neutron flux is
j HrL =
S
4 p D
r
where D is the diffusion constant, L = D ê S a is called the diffusion length (refer to the class notes). Develop a Monte Carlo
simulation for this system and show that the flux follows this form, evaluate the diffusion length, and compare with the analytic
result. (One strategy is to launch a neutron from the origin in an arbitrary direction and track the trajectory until it is absorbed.
Assume that the neutron scatters isotropically from the nuclei of the moderating material. Keep track of the radius of the neutron
(suitably binned) for the purposes of constructing j, calculating the effective diffusion length at the end of the simulation, and
compare your Monte Carlo result with the analytic result, Eq.(7).)
ü Analytic results for mean free path and diffusion lengths
H* Use Sodium thermal neutron data *L
Avalue = 23 ;
data = 8 r Ø .00097, n d Ø .02541, Ss Ø .08131, Sg Ø .01347 , Sf Ø 0 , n Ø 0 < ;
H* Use Carbon thermal neutron data *L
Avalue = 12 ;
C12data = 8 r Ø .00160, n d Ø .08023, Ss -> .3811, Sg Ø .0002728 , Sf Ø 0 , n Ø 0 < ;
muavg @ A _D = 2 ê H 3 A L ;
Sa = HSg + SfL ê. data;
MFP = 1 ê Ss ê. data ê. Cons;
MFPabs = 1 ê Sa ê. data ê. Cons;
MFPtotal = 1 ê HSs + Sg + SfL ê. data ê. Cons;
Diff = 1 ê HSa + Ss H 1 - muavg @ Avalue DLL ê 3 ê. data ê. Cons;
DAlt = Ss ê HSa + Ss H 1 - muavg @ Avalue DLL ^ 2 ê 3 ê. data ê. Cons;
L = Sqrt @ Diff ê SaD ê. data ê. Cons;
LAlt = Sqrt @ DAlt ê SaD ê. data ê. Cons;
AbsorptionLength = 1 ê Sa ê. data ê. Cons;
DiffTh = vspeed Diff;
H Ss* , Sa Ø cm^ - 1 *L
Print @ " Mean Free Path H scattering L = ", MFP, " cm" D
Print @ " Mean Free Path H absorption L = ", MFPabs, " cm" D
Print @ " Mean Free Path H total L = ", MFPtotal, " cm" D
Print @ " Flux Diffusion constant: D = ", Diff, " cm,", " D H alt L = ", DAlt, " cm" D
Print @ " Diffusion Length: L = ", L, " cm,", " L H alt L = ", LAlt, " cm" D
Print @ " Absorption MFP: = ", AbsorptionLength, " cm" D
Print @ " Theoretical Density Diffusion constant = ", DiffTh, " cm^ 2 ê s" D
Mean Free Path HscatteringL = 12.2986 cm
Mean Free Path HabsorptionL = 74.239 cm
Mean Free Path HtotalL = 10.5507 cm
Flux Diffusion constant: D = 3.69292 cm, DHaltL = 3.32664 cm
Diffusion Length: L = 16.5577 cm, LHaltL = 15.7152 cm
Absorption MFP: = 74.239 cm
Theoretical Density Diffusion constant = 916 073. cm^ 2 ês
ü Monte Carlo tracking of neutron paths
ü Monte Carlo simulation
Nexp = number of "experiments"
Nneutrons= number of neutrons per experiment
dt = 5. 10 ^ - 6 ; ds = vspeed dt;
dPs = Ss ds ê. data ê. Cons;
dPa = HSg + SfL ds ê. data ê. Cons;
Print @ " For time step of ", dt, " sec, and average speed of, v = ",
vspeed, " cm ê s, distance per step = ", ds, " cm" D
Print @ " The probability of scattering in one time step:" D
Print @ " Pscatt = ", dPs, " Pabs = ", dPa D
For time step of 5. μ 10
sec, and average speed of, v =
248 062. cmês, distance per step = 1.24031 cm
The probability of scattering in one time step:
Pscatt = 0.10085 Pabs = 0.
Clear @ LenCos D
Nexp = 100 ; Nneutrons = 1000 ;
Rmax = 4 L; NRbins = 40 ; dR =
Rmax
NRbins
; NumDen = Table @ 0 , 8 i, 1 , NRbins + 1 <D ;
vcmvalue = vcm @ Avalue D ; voverAplus1 =
vspeed
1 + Avalue
Timing B CosAvgTable = 8 < ; MFPAvgTable = 8 < ; RAvgTable = 8 < ; R2AvgTable = 8 < ; RAbsAvgTable = 8 < ;
NAbsAvgTable = 8 < ; Do B Rtable = 8 < ; MFPtable = 8 < ; costable = 8 < ; Rabstable = 8 < ; nabstable = 8 < ;
Do B For B 8 vhat = 8 0 , 0 , 1 < ; r = 8 0 , 0 , 0 < ; NoScatt = 0 ; NoScattList = 8 < ; cossum = 0 ; nscatt = 0 ; iStep = 1 ; iStop = - 1 < ,
iStop < 0 && iStep < 10
5 , iStep ++ , : ran = RandomReal @D ;
If B ran > dPs + dPa, NoScatt ++ , If B ran < dPs, : vhatold = vhat; thcm = p RandomReal @D ; phicm = 2 p RandomReal @D ;
vlab = vcmvalue 8 Sin @ thcm D Cos @ phicm D , Sin @ thcm D Sin @ phicm D , Cos @ thcm D< + voverAplus1 vhat;
vhat =
vlab
vlab.vlab
; nscatt ++ ; cos = vhat.vhatold; cossum = cossum + cos;
AppendTo @ NoScattList, NoScatt D ; NoScatt = 0 > , AppendTo @ NoScattList, NoScatt D ;
AppendTo @ nabstable, iStep D ; AppendTo @ Rabstable, rmag D ; iStop = 1 FF ; r = r + ds vhat; rmag = r.r ;
rindex = Floor B
rmag
dR
F + 1 ; If @ rindex > NRbins, NumDen P NRbins + 1 T ++ , NumDen P rindex T ++D ; >F ;
If B nscatt > 0 , AppendTo B costable,
cossum
nscatt
FF ; AppendTo B Rtable, r.r F ; LenNos = Length @ NoScattList D ;
If B LenNos > 0 , AppendTo B MFPtable, N B
ds ⁄i= 1
LenNos NoScattList P i T
LenNos
FFF ;, 8 j, 1 , Nneutrons <F ;
LenCos = Length @ costable D ; AppendTo B CosAvgTable, N B
⁄i= 1
LenCos costable P i T
LenCos
FF ;
LenMFP = Length @ MFPtable D ;
AppendTo B MFPAvgTable, N B
⁄i= 1
LenMFP MFPtable P i T
LenMFP
FF ;
LenR = Length @ Rtable D ; AppendTo B RAvgTable, N B
⁄i= 1
LenR Rtable P i T
LenR
FF ;
AppendTo B R2AvgTable, N B
⁄i= 1
LenR Rtable P i T
2
LenR
FF ;
LenRabs = Length @ Rabstable D ;
AppendTo B RAbsAvgTable,
⁄i= 1
LenRabs Rabstable P i T
LenRabs
F ; LenNabs = Length @ nabstable D ;
AppendTo B NAbsAvgTable,
⁄i= 1
LenNabs nabstable P i T
LenNabs
F ;, 8 iexp, 1 , Nexp <F ; F
8 287.685, Null<
ü Plot results
PhiTable = Table B If B NumDen P i T ã 0 ,
Rmax
2
NumDen P i T
H i dR L
2
F , 8 i, 1 , NRbins <F ;
RPhiTable = Table @ i dR PhiTable P i T , 8 i, 3 , NRbins <D ;
LogRPhiTable = Table @ Log @ RPhiTable P i TD , 8 i, 1 , Length @ RPhiTable D<D ;
PhiFit @ x _D = Fit @ LogRPhiTable, 8 1 , x < , x D ;
LFit = -
dR
PhiFit @ 1 D - PhiFit @ 0 D
; LMC =
Ravg
Print @ " Diffusion Length Monte Carlo L = ", LMC, " cm" D
Print @ " Fit to simulation data L = ", LFit, " cm" D
Print @ " Diffusion L = ", L, " cm" D
Print @ " Diffusion - alt L = ", LAlt, " cm" D
horizaxis = Style @ "r H dR units L ", FontFamily Ø "Tahoma", FontColor Ø Blue, FontWeight Ø Bold, FontSize Ø 12 D ;
vertaxis = Style @ "Log @ r jD ", FontFamily Ø "Tahoma", FontColor Ø Blue, FontWeight Ø Bold, FontSize Ø 12 D ;
plotname = Style @ "Fit to Log @ r jD ", FontFamily Ø "Tahoma", FontColor Ø Black, FontWeight Ø Bold, FontSize Ø 14 D ;
p1 = ListPlot @ LogRPhiTable, DisplayFunction Ø Identity,
PlotStyle Ø 8 RGBColor @ 1 , 0 , 0 D , PointSize @ 0.02` D< , Frame Ø True, GridLines Ø Automatic, PlotLabel Ø plotname,
FrameLabel Ø 8 horizaxis, vertaxis < , ImageSize Ø 400 , Background Ø LightOrange, PlotRange Ø All D ;
p2 = Plot @ PhiFit @ x D , 8 x, 0 , NRbins < , DisplayFunction Ø Identity, PlotStyle Ø 8 Blue, Thickness @ 0.005` D< ,
Frame Ø True, GridLines Ø Automatic, PlotLabel Ø plotname, FrameLabel Ø 8 horizaxis, vertaxis < ,
ImageSize Ø 400 , Background Ø LightOrange, PlotRange Ø All D ;
Show @ p1, p2, DisplayFunction Ø $DisplayFunction D
Diffusion Length Monte Carlo L = 15.3992 cm
Fit to simulation data L = 15.5973 cm
Diffusion L = 16.5577 cm
Diffusion-alt L = 15.7152 cm
r HdR unitsL
Log
@
r
j
D
Fit to Log@r jD
norm = (^) ‚
i= 1
NRbins
NumDen P i T ; DenNormed = Table B N B
NumDen P i T
norm
F , 8 i, 1 , NRbins <F ;
ShellVol @ iR _D = dR
3 4 p iR
2 ; PhiMCTable = Table B: i -
dR,
DenNormed P i T
ShellVol @ i D
> , 8 i, 1 , NRbins <F ;
Phi @ x _D =
-
x LD
x
; norm = 4 p (^) ‡ 0
¶
Phi @ x D x
2 „ x; Phi0 @ x _D =
Phi @ x D
norm
ê. LD Ø L;
PhiAnalyticTable = Table B: i -
dR, Phi0 B i -
dR F> , 8 i, 1 , NRbins <F ;
horizaxis = Style @ "r H cm L ", FontFamily Ø "Tahoma", FontColor Ø Blue, FontWeight Ø Bold, FontSize Ø 12 D ;
vertaxis = Style @ " j ", FontFamily Ø "Tahoma", FontColor Ø Blue, FontWeight Ø Bold, FontSize Ø 12 D ;
plotname = Style @ "Neutron Flux", FontFamily Ø "Tahoma", FontColor Ø Black, FontWeight Ø Bold, FontSize Ø 14 D ;
p1 = ListPlot @ PhiMCTable, DisplayFunction Ø Identity,
PlotStyle Ø 8 RGBColor @ 1 , 0 , 0 D , PointSize @ 0.02` D< , Frame Ø True, GridLines Ø Automatic, PlotLabel Ø plotname,
FrameLabel Ø 8 horizaxis, vertaxis < , ImageSize Ø 400 , Background Ø LightOrange, PlotRange Ø All D ;
p2 = ListPlot @ PhiAnalyticTable, Joined Ø True, DisplayFunction Ø Identity,
PlotStyle Ø 8 Blue, Thickness @ 0.004` D< , Frame Ø True, GridLines Ø Automatic, PlotLabel Ø plotname,
FrameLabel Ø 8 horizaxis, vertaxis < , ImageSize Ø 400 , Background Ø LightOrange, PlotRange Ø All D ;
Show @ p1, p2, DisplayFunction Ø $DisplayFunction D
Task 2: Criticality for the infinite reactor
Task 2. Consider an infinite homogeneous mix of fissile
U (10% by weight) and sodium (90% by weight) at room tempera-
ture. (a) First find the composition average values for the macroscopic scattering, absorption, and fission cross sections. (b)
Develop a Monte Carlo simulation for this reactor to estimate the neutron multiplication factor, k¶, and compare with the
expected analytic result,
k¶ =
n Sf
Sa
Track the neutron locations and plot the Monte Carlo modeled flux as a function of radius (suitably binned and normalized).
ü Calculate the macroscopic cross sections for the mixture (10% 235U, 99% Na) (by weight) and
analytic value of k¶.
Clear @ soln, massNa, nNa, massU, nU, VNa, VU, Mix D
H* Thermal neutron values *L
Nadata = 8 r Ø .00097, n d Ø .02541, Ss Ø .08131, Sg Ø .01347 , Sf Ø 0 , n Ø 0 < ;
U235data = 8 r Ø .01886, n d Ø .04833, Ss Ø .01588, Sg Ø 4.833, Sf -> 28.37, n Ø 2.42 < ;
VNa = 1 ë n d ê. Nadata; VU = 1 ë n d ê. U235data;
amu = 1.67 μ 10 ^ - 27 ; rhoNa = r ê. Nadata; rhoU = r ê. U235data;
Mix = .1;
soln = Flatten @ Solve @ 8 nNa massNa Mix ã H 1 - Mix L massU nU, nNa VNa + nU VU == 1 < , 8 nNa, nU <DD ;
NdenNa = HH nNa ê. soln L ê. 8 VNa Ø massNa ê rhoNa, VU Ø massU ê rhoU <L ê. 8 massU Ø 235 amu, massNa Ø 23 amu < ;
NdenU = HH nU ê. soln L ê. 8 VNa Ø massNa ê rhoNa, VU Ø massU ê rhoU <L ê. 8 massU Ø 235 amu, massNa Ø 23 amu < ;
H* check *L
NdenNa massNa ê H NdenU massU L ê. 8 massU Ø 235 amu, massNa Ø 23 amu < ;
Print @ " Number density of Sodium = ", NdenNa, " H 10 ^ 24 ê cm^ 3 L " D
Print @ " Number density of Uranium = ", NdenU , " H 10 ^ 24 ê cm^ 3 L " D
SigGMix = II NdenNa ë n dM Sg ê. Nadata M + II NdenU ë n dM Sg ê. U235data M ;
SigSMix = II NdenNa ë n dM Ss ê. Nadata M + II NdenU ë n dM Ss ê. U235data M ;
SigFMix = II NdenNa ë n dM Sf ê. Nadata M + II NdenU ë n dM Sf ê. U235data M ;
NdenMix = NdenNa + NdenU;
rhoMix = H NdenNa 23 amu + NdenU 235 amu L ;
etaMix = n ê. U235data;
data = 8 r Ø rhoMix, n d Ø NdenMix, Ss Ø SigSMix, Sg Ø SigGMix, Sf Ø SigFMix, n Ø etaMix <
kinf = Hn Sf ê HSf + SgLL ê. data;
Print @ " k ¶ H Analytic L = ", kinf D
Timing B Nneutrons = 20 000 ; Rtable = 8 < ; vcmvalue = vcm @ Avalue D ; voverAplus1 =
v
1 + Avalue
nfiss = 0 ; Do B For B 8 vhat = 8 0 , 0 , 1 < ; r = 8 0 , 0 , 0 < ; NoScatt = 0 ; iStep = 1 ; iStop = - 1 < ,
iStop < 0 && iStep < 10
4 , iStep ++ , : ran = RandomReal @D ; If B ran > dPs + dPg + dPf, Null,
If B ran < dPs, : vhatold = vhat; thcm = p RandomReal @D ; phicm = 2 p RandomReal @D ; vlab =
vcmvalue 8 Sin @ thcm D Cos @ phicm D , Sin @ thcm D Sin @ phicm D , Cos @ thcm D< + voverAplus1 vhat; vhat =
vlab
vlab.vlab
If @ ran < dPs + dPf, nfiss ++D ; iStop = 1 FF ; r = r + ds vhat; >F ; AppendTo B Rtable, r.r F ;, 8 j, 1 , Nneutrons <F ; F
8 16.0877 Second, Null<
kinfMC = Hn nfiss ê Nneutrons L ê. data;
Print @ " k ¶ H Monte Carlo L = ", kinfMC, " H Analytic L = ", kinf D ;
H* Average radius *L
LenR = Length @ Rtable D ;
Ravg = N @ Sum @ Rtable @@ i DD , 8 i, 1 , LenR <D ê LenR D ;
Print A " Diffusion length H MC = R
ê ê 2 L = ", Ravg ê 2 , " cm Analytic = ", L, " cm" E
k¶ HMonte CarloL = 1.93116 HAnalyticL = 1.
Diffusion length HMC = R
ê
ê 2 L = 2.15267 cm Analytic = 2.42351 cm
