Time-dependent Monte Carlo in fissile systems with beta-delayed neutron precursors
TTime-dependent Monte Carlo in fissilesystems with beta-delayed neutronprecursors
Tesisentregada a laUniversidad de Chileen cumplimiento parcial de los requisitospara optar al grado deDoctor en Ciencias con menci´on en F´ısicaFacultad de Ciencias F´ısicas y Matem´aticaspor
Jaime Alfonso Romero Barrientos
Enero, 2021Director de Tesis:
Dr. Hugo Arellano Sep´ulveda
Co-Director de Tesis:
Dr. Francisco Gabriel Molina Palacios a r X i v : . [ phy s i c s . c o m p - ph ] J a n LOSSARY β -delayed neutron emitter Nuclei that emits a β -delayed neutron. Analog Monte Carlo simulation
Monte Carlo simulation which followsthe natural (i.e. physical) probability distribution function for randomsampling (See Sec. 2.4.3).
Batch
In a Monte Carlo simulation, a batch is a single realization of a tallyrandom variable. In the simulation both the number of batches as well asthe number of particles per batch must be specified.
Configuration
In this work, it refers to the state of neutron multiplicationof the system. Configuration can be either subcritical, critical or supercrit-ical.
Criticality calculation
Also called eigenvalue calculation. Monte Carlosimulation where the objective is the determination of the state of neutronmultiplication in a fissile system. If k eff = 1 the system is in a criticalconfiguration. While if k eff < k eff >
1) the system is in a subcritical(supercritical) configuration (See Sec. 2.1.3).
Delayed neutron fraction
Denoted as β . Represents the fraction of to-tal fission neutrons which are delayed (See Sec. 2.2). Effective delayed neutron fraction
Denoted as β eff . Represents thefraction of fissions caused by delayed neutrons. Its defined as the delayedneutron fraction β weighted by the neutron importance, which representshow effective the neutron is in causing fission (See Sec. 2.1.4). Effective multiplication factor
Parameter that shows the state of neu-tron multiplication in a fissile system. If k eff = 1 the system is in a criticalconfiguration. While if k eff < k eff >
1) the system is in a subcritical(supercritical) configuration (See Sec. 2.1.3).
Fixed-source calculation
Monte Carlo simulation where the initial par-ticle source is known and the resulting neutron distribution is determined.
CNP
MCNP is a general-purpose Monte Carlo N-Particle code that canbe used for neutron, photon, electron, or coupled neutron/photon/electrontransport. Developed and mantained by Los Alamos National Laboratory. N -group structure Scheme that groups all the β -delayed neutron precur-sors into N -groups. Each precursor group contains a number of differentisotopes. In ENDF/B.VIII.0 database, N = 6. In JEFF-3.1.1, N = 8 (SeeSec. 2.2). Non-analog Monte Carlo simulation
Monte Carlo simulation which fol-lows a modified (i.e. non-physical) probability distribution function for ran-dom sampling in order to reduce the variance of the result obtained whenusing an analog Monte Carlo simulation (See Sec. 2.4.3).
OpenMC
OpenMC is a community-developed Monte Carlo neutron andphoton transport simulation code. It is capable of performing fixed source,k-eigenvalue, and subcritical multiplication calculations on models builtusing either a constructive solid geometry or CAD representation. It wasoriginally developed by members of the Computational Reactor PhysicsGroup at the Massachusetts Institute of Technology starting in 2011 (SeeSec. 3.1).
OpenMC(TD)
Time-dependent OpenMC. OpenMC code with added ca-pabilities shown in this work, that is: explicit presence of time, scoring oftime dependent quantities, non-analog simulation scheme to simulate the β -delayed neutron emission, population control to keep the number of parti-cles constant and option to use either the N -group precursor structure (with N = 1 , , or 8) or M -individual precursor structure (with 1 ≤ M ≤ pcm Per-cent mille one-thousandth of a percent.
Point kinetics approximation
Theoretical approximation used to studythe kinetic behaviour of a fissile system, where the flux is assumed to be aseparable function of space and time (See Sec. 2.1.4).
Prompt drop
Fast decrease in neutron population or flux caused by aeduction in the system reactivity. The timescale of this process is of theorder of the prompt neutron generation time (Sec. 4.3.1).
Prompt jump
Fast increase in neutron population or flux caused by anincrease in the system reactivity. The timescale of this process is of the orderof the prompt neutron generation time (See Sec. 4.2.3 and Sec. 4.3.2).
Prompt neutron generation time At k eff = 1, average time betweentwo generations of prompt neutrons. Precursor structure
In this work, it refers to the organization of the β -delayed neutron emitters. Precursor structure can be either a N -groupprecursor structure or a N -individual precursor structure. In the latter, β -delayed neutrons are emitted from individual precursors, not groups (SeeSec. 3.3.2 and Sec. 4.4). Precursor
Fission product (
Z, N ) which decays through a β -delayed pro-cess to another nuclei ( Z +1 , N − Z +1 , N − β -delayed neutron (See Sec. 2.2). Skipped cycles
Batches discarded before scores begin to accumulate in aMonte Carlo calculation.
System
In this work, it refers to the simulated structure, characterized byits geometry, materials, moderation, and so on.
Transient source
In this work, it refers to the initial particle source, com-prised of neutrons and precursors, needed to start a transient calculation(See Sec. 3.4).
Weight, statistical
Number that represents how many real (i.e. physical)particles a Monte Carlo particle represents. If the statistical weight of aneutron is 2, then that neutron represents 2 neutrons.
ESUMEN
En el campo de la f´ısica de reactores nucleares, los fen´omenos transientessuelen estudiarse usando m´etodos deterministas o h´ıbridos. Estos m´etodos requierende variadas aproximaciones, tales como: discretizaciones de la geometr´ıa, del tiempoy de la energ´ıa; homogeneizaci´on de materiales; y suposici´on de condiciones de di-fusi´on, por mencionar algunas. En este contexto, las simulaciones Monte Carlo sonespecialmente adecuadas para estudiar estos problemas. Los retos que se presentanal usar simulaciones Monte Carlo en cin´etica espacio-temporal de sistemas fisibles sonlas escalas de tiempo inmensamente distintas involucradas en la emisi´on de neutronesinmediatos y retardados, lo que implica que los resultados obtenidos tienen asociadauna gran varianza si se utiliza una simulaci´on Monte Carlo an´aloga. Adem´as, tantoen simulaciones deterministas como en Monte Carlo, los precursores de neutronesretardados est´an agrupados en una estructura 6 u 8 grupos, pero hoy en d´ıa no hayuna raz´on s´olida para mantener esta agrupaci´on.En este trabajo, y por primera vez, se han implementado los datos deprecursores individuales en una simulaci´on Monte Carlo, incluyendo expl´ıctamente ladependencia temporal relacionada con la emisi´on β retardada de neutrones. Esto fuelogrado modificando el c´odigo abierto Monte Carlo OpenMC. En el c´odigo modificado– Time Dependent
OpenMC u OpenMC(TD) – se abord´o la dependencia temporalrelacionada con la emisi´on retardada de neutrones originada de la desintegraci´on β .La varianza del valor esperado de observables, como el flujo neutr´onico, asociada a lasdiferentes escalas de tiempo entre los neutrones inmediatos y prompts, fue reducidaforzando la desintegraci´on de una nueva part´ıcula Monte Carlo a˜nadida al c´odigo,el precursor, dentro de cada intervalo temporal, incrementando intencionalmente eln´umero de neutrones retardados en la simulaci´on. Dado que hay una producci´oncontinua de neutrones retardados, se tuvo que imponer el control de poblaci´on. Estose logr´o usando el m´etodo de combing al final de cada intervalo temporal.Los datos de secciones eficaces dependientes de la energ´ıa vienen de labiblioteca JEFF-3.1.1. Los datos de los precursores individuales fueron tomados delas bibliotecas JEFF-3.1.1 ( yields cumulativos) y ENDF-B/VIII.0 (probabilidades deemisi´on de neutrones retardados y espectros de energ´ıa de neutrones retardados).penMC(TD) fue probado en: i) un sistema monoenerg´etico; ii) un sis-tema sin moderaci´on y dependiente de la energ´ıa donde los precursores se tomaronindividualmente o en grupos; y finalmente iii) un sistema moderado por agua li-viana dependiente de la energ´ıa, usando 6-grupos de precursores, 50 precursores y40 precursores individuales. ABSTRACT
In the field of nuclear reactor physics, transient phenomena are usuallystudied using deterministic or hybrids methods. These methods require many ap-proximations, such as: geometry, time and energy discretizations, material homog-enization and assumption of diffusion conditions, among others. In this context,Monte Carlo simulations are specially adequate to study these problems. Challengespresented when using Monte Carlo simulations in space-time kinetics in fissile sys-tems are the immensely different time-scales involved in prompt and delayed neutronemission, which implies that results obtained have a large variance associated if ananalog Monte Carlo simulation is utilized. Furthermore, in both deterministic andMonte Carlo simulations delayed neutron precursors are grouped in a 6- or 8- groupstructure, but nowadays there is not a solid reason to keep this aggregation.In this work, and for the first time, individual precursor data is implementedin a Monte Carlo simulation, explicitly including the time dependence related to the β -delayed neutron emission. This was accomplished by modifying the open sourceMonte Carlo code OpenMC. In the modified code – Time Dependent OpenMC orOpenMC(TD) – time dependency related to delayed neutron emission originatedfrom β -decay was addressed. The variance of the expected values of observables, suchas neutron flux, associated to the different time scales between prompt and delayedneutrons was reduced by forcing the decay of a new Monte Carlo particle-like addedto the code, the precursor, within each time interval, intentionally increasing thenumber of delayed neutrons in the simulation. Since there is a continuous productionof delayed neutrons, population control had to be enforced. This was accomplishedby using the combing method at the end of each time interval.Continuous energy neutron cross-sections data used comes from JEFF-3.1.1library. Individual precursor data was taken from JEFF-3.1.1 (cumulative yields) andENDF-B/VIII.0 (delayed neutron emission probabilities and delayed neutron energyspectra). OpenMC(TD) was tested in: i) a monoenergetic system; ii) an energy de-pendent unmoderated system where the precursors were taken individually or in agroup structure; and finally iii) a light-water moderated energy dependent system,using 6-groups, 50 and 40 individual precursors. ontents U system . . . . . . . . . . . . . . . . . . . . . . 514.3.1 Subcritical configuration . . . . . . . . . . . . . . . . . . . . . 534.3.2 Supercritical configuration . . . . . . . . . . . . . . . . . . . . 634.4 Light-water moderated energy dependent system with individual pre-cursor structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 764.4.1 Criticality calculation using individual precursors . . . . . . . 784.4.2 Critical configuration with 50 individual precursor structure . 794.4.3 Critical configuration without the 10 most important precursors 82 i C Monoenergetic fissile system with 1-group precursor structure 97D Summary of the simulations performed in this work 99E Individual precursor data 102 ist of Figures
U. Data was retrieved from the JEFF-3 . . Br delayed neutron precursor . . . . . . . . . . . . . . 122.3 Schematic representation of the prompt and β -delayed neutron emission. 133.1 Schematic representation of the time scales associated to the delayedneutron emission and the lifetime of the prompt chains. This differenttime scales produce large variance in the quantities scored. . . . . . . 273.2 Schematic representation of the forced decay scheme for the precur-sors, where the precursor is forced to decay at the beginning of eachtime bin, so there are scored in every time interval. . . . . . . . . . . 333.3 Diagram of the application of the combing method for 4 particles oftotal weight W combed into M = 3. The particles kept by the combare particle 1, particle 3 and particle 4, each with weight W/
3. . . . . 354.1 Time evolution of the neutron flux in a subcritical configuration forthe RA-6 reactor, obtained from running a fixed source simulation inboth OpenMC(TD) and MCNP. Results were scored for 10 ms. Bothcodes are in good agreement with the experimental benchmark result. 39iiiv4.2 Neutron flux as a function of time in a simple transport problem. Inred, the neutron flux obtained for a non-transient fixed source simu-lation is shown. In blue the neutron flux obtained from a transientsimulation divided in time intervals is shown. Both results are equiv-alent. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.3 Time evolution of the neutron flux for the monoenergetic system stud-ied. Results obtained using OpenMC(TD) are shown in blue, whilethe fit to the point kinetics solution given by Eq. (4.2) is shown in red. 424.4 β -delayed neutron activity for 6 precursor groups and 50 individualprecursors is shown. In blue, A ( t ) is shown, while A ( t ) is shown inred. As it can be seen, both activities are equivalent. . . . . . . . . . 444.5 Time evolution of neutron flux for a subcritical configuration with k eff = 0 . ± . µ s each.Results of neutron flux when delayed neutron energy was sampledfrom spectra are shown in red (behind the blue curve), while resultsobtained when using the average energy for delayed neutron emissionare shown in blue. It can be seen that both results are equivalent. . . 464.6 Time evolution of the neutron population for a monoenergetic systemin a subcritical configuration ( k eff = 0 . ± . k eff =1 . ± . t = 10 s a reactivity of 211 pcm is inserted. After30 s the system is brought back to critical configuration. . . . . . . . 504.9 Group 1 β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 534.10 Study i). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, and then the decay of the neutron populationslows. Inset figure shows the prompt drop zoomed for the first 5 µ s. . 544.11 Study ii). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, and then the decay of the neutron populationslows. Inset figure shows the prompt drop zoomed for the first 5 µ s. . 564.12 Study iii). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, and then the decay of the neutron populationslows. Inset figure shows the prompt drop zoomed for the first 5 µ s. . 58i4.13 Study iv). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, and then the decay of the neutron populationslows. Inset figure shows the prompt drop zoomed for the first 5 µ s. . 594.14 Study v). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, and then the decay of the neutron populationslows. Inset figure shows the prompt drop zoomed for the first 5 µ s. . 614.15 Study vi). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, and then the decay of the neutron populationslows. Inset figure shows the prompt drop zoomed for the first 5 µ s. . 62ii4.16 Study i). Time evolution of the neutron flux in the studied super-critical configuration. The initial transient source is prepared in acritical state and at the beginning of the transient Monte Carlo simu-lation using OpenMC(TD) the configuration is made supercritical byincreasing δ U . The time evolution of the neutron flux is shown inblue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutronpopulation slows. Inset figure shows the prompt jump zoomed for thefirst 10 µ s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 654.17 Study ii). Time evolution of the neutron flux in the studied super-critical configuration. The initial transient source is prepared in acritical state and at the beginning of the transient Monte Carlo simu-lation using OpenMC(TD) the configuration is made supercritical byincreasing δ U . The time evolution of the neutron flux is shown inblue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutronpopulation slows. Inset figure shows the prompt jump zoomed for thefirst 10 µ s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 664.18 Study iii). Time evolution of the neutron flux in the studied super-critical configuration. The initial transient source is prepared in acritical state and at the beginning of the transient Monte Carlo simu-lation using OpenMC(TD) the configuration is made supercritical byincreasing δ U . The time evolution of the neutron flux is shown inblue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutronpopulation slows. Inset figure shows the prompt jump zoomed for thefirst 10 µ s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68iii4.19 Study iv). Time evolution of the neutron flux in the studied super-critical configuration. The initial transient source is prepared in acritical state and at the beginning of the transient Monte Carlo simu-lation using OpenMC(TD) the configuration is made supercritical byincreasing δ U . The time evolution of the neutron flux is shown inblue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutronpopulation slows. Inset figure shows the prompt jump zoomed for thefirst 10 µ s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 704.20 Study v). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, and then the decay of the neutron populationslows. Inset figure shows the prompt drop zoomed for the first 10 µ s. 724.21 Study vi). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a criticalstate and at the beginning of the transient Monte Carlo simulationusing OpenMC(TD), the system is made subcritical by decreasing δ U . The time evolution of the neutron flux is shown in blue, whilethe fit obtained is shown in red. Between 0 < t < µ s the promptdrop can be observed, after which the decay of the neutron flux slows.Inset figure shows the prompt drop zoomed for the first 10 µ s. . . . . 744.22 Time evolution of neutron flux in a water moderated box made ofpure U, simulated with OpenMC. Results for 6 groups are shownin blue, and results for 50 individual precursors are shown in red.Both results show a slightly supercritical system, where neutron fluxincreases slowly in time, consistent with the k eff of a near criticalconfiguration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80x4.23 Time evolution for the neutron flux in a water moderated box madeof pure U, simulated with OpenMC. When the 6 groups are usedresults shown in blue. Results obtained when 50 individual precursorsare used are shown in red, and results obtained when using 40 pre-cursors are shown in green. In this case it can be seen that the timeevolution of the neutron flux calculated using 40 precursors clearly di-verges from the previous results, showing that the neutron flux growsmore rapidly in time. . . . . . . . . . . . . . . . . . . . . . . . . . . . 83A.1 Group 1, β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 90A.2 Group 2, β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 91A.3 Group 3, β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 91A.4 Group 4, β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 92A.5 Group 5, β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 92A.6 Group 6, β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 93A.7 Group 7 β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 93A.8 Group 8 β -delayed neutron energy spectrum from JEFF-3.1.1. . . . . 94C.1 Box of 10 cm ×
12 cm ×
20 cm simulated in Sec. 4.2 and Sec. 4.3. . . 97 ist of Tables
U fission. In 6-th groupseveral isotopes of Br, As, Rb are included. . . . . . . . . . . . . . . 122.2 Independent yields, cumulative yields and branching ratio values foundin JEFF 3 . . . β -delayed neutronemission using the precursor group structure or the individual pre-cursors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304.1 Decay constants obtained for the time evolution of the neutron fluxobtained using the pulsed method in the RA-6 reactor. . . . . . . . . 404.2 Decay constants obtained for the time evolution of the neutron fluxusing the RA-6 reactor. . . . . . . . . . . . . . . . . . . . . . . . . . . 424.3 Results obtained for the reactivity of the monoenergetic simulatedsystem in a subcritical configuration using 1-group precursor structure. 48xi4.4 Results obtained for the reactivity of the monoenergetic simulatedsystem in a critical configuration using 1-group precursor structure. . 504.5 Study i). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a subcritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting Eq. (B.4) to thetime evolution of the neutron flux. . . . . . . . . . . . . . . . . . . . . 554.6 Study ii). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a subcritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting Eq. (B.4) to thetime evolution of the neutron flux. . . . . . . . . . . . . . . . . . . . . 564.7 Study iii). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a subcritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting Eq. (B.4) to thetime evolution of the neutron flux. . . . . . . . . . . . . . . . . . . . . 574.8 Study iv). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a subcritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting Eq. (B.4) to thetime evolution of the neutron flux. . . . . . . . . . . . . . . . . . . . . 604.9 Study v). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a subcritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting Eq. (B.4) to thetime evolution of the neutron flux. . . . . . . . . . . . . . . . . . . . . 61ii4.10 Study vi). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a subcritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting Eq. (B.4) to thetime evolution of the neutron flux. . . . . . . . . . . . . . . . . . . . . 634.11 Study i). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a supercritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting point kinetics so-lution to the time evolution of the neutron flux. . . . . . . . . . . . . 644.12 Study ii). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a supercritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting point kinetics so-lution to the time evolution of the neutron flux. . . . . . . . . . . . . 674.13 Study iii). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a supercritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting point kinetics so-lution to the time evolution of the neutron flux. . . . . . . . . . . . . 694.14 Study iv). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a supercritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting point kinetics so-lution to the time evolution of the neutron flux. . . . . . . . . . . . . 714.15 Study v). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a supercritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while theOpenMC(TD) parameters were obtained by fitting point kinetics so-lution to the time evolution of the neutron flux. . . . . . . . . . . . . 73iii4.16 Study vi). Values of the parameters obtained from running an OpenMC(TD)transient simulation in a supercritical configuration. The calculatedvalues for the parameters were obtained using MCNP, while the OpenMC(TD)parameters were obtained by fitting point kinetics solution to the timeevolution of the neutron flux. . . . . . . . . . . . . . . . . . . . . . . 754.17 Effective multiplication factors obtained for the U cube when isthermalized by surrounding it with a water moderator of a 4 .
29 cmthickness. Results for 6-group structure and 50 individual precursors. 784.18 Results obtained for the reactivity of the water moderated energydependent simulated system in a critical configuration using 6-groupand 50 individual precursor structure. . . . . . . . . . . . . . . . . . . 814.19 Results obtained for the reactivity of the water moderated energydependent simulated system in a critical configuration using 6-group,50 individual and 40 individual precursor structures. . . . . . . . . . 834.20 Summary of all results obtained with OpenMC(TD).(*) Error wasobtained only from the calculated result using point kinetics equationsor criticality calculation using OpenMC, given that OpenMC(TD)error could be improved as explained in Sec. 4.4.2. . . . . . . . . . . . 84C.1 Material cross sections and parameters of the monoenergetic system. . 98D.1 Summary of simulation parameters for monoenergetic fissile systemin subcritical, critical, and reactivity insertion configurations, using1-group precursor structure. . . . . . . . . . . . . . . . . . . . . . . . 100D.2 Summary of simulation parameters for energy dependent system in asubcritical configuration, using different precursor structures. . . . . . 100D.3 Summary of simulation parameters for energy dependent system in asupercritical configuration, using different precursor structures. . . . . 101ivD.4 Summary of simulation parameters for light-water moderated, energydependent system in a critical configuration, using 6-group, 40 indi-vidual, and 50 individual precursor structures. . . . . . . . . . . . . . 101E.1 6-group precursor structure used in this work. . . . . . . . . . . . . . 102E.2 50 individual precursor structure used in this work. Precursors areordered by importance. . . . . . . . . . . . . . . . . . . . . . . . . . . 104E.3 40 individual precursor structure used in this work. Precursors areordered by importance. . . . . . . . . . . . . . . . . . . . . . . . . . . 105 hapter 1Introduction
Nuclear fission is a process where the atomic nucleus splits in two or three fissionproducts (lighter weight nuclei) and neutrons. In heavy nuclei this process canhappen as a spontaneous desintegration (
Cf), or it can be induced by the reactionwith a neutron. If fission is induced in a nucleus by a thermal energy neutron, thenthe nucleus is said to be fissile (
U or
Pu). If the nucleus requires neutrons witha certain threshold energy to be fissioned, then it is said to be a fisssionable nucleus(
U or
Th).In fission reactions two types of neutrons are emitted: prompt and delayedneutrons. Prompt neutrons are emitted almost instantaneously ( ∼ − s) afterfission occurs with energies of the order of a few MeV. Meanwhile, delayed neutronsare emitted from milliseconds to tens of seconds after fission with energies of theorder of hundreds of keV. Delayed neutron emission is associated to the decay ofisotopes from the decay chain of fission products. These nuclei emitters of β -delayedneutrons, are called precursors. For example, for U there are about 540 fissionproducts and 270 precursors [1].If there is enough fissile material, a neutron can induce fission in other nu-cleus and initiate a chain reaction. This chain reaction can be sustained in timedepending on the density of fissile material, neutron energy at the moment of fissionand the fission reaction rate. The Neutron Transport Equation models the propaga-tion of neutrons in a fissile system. 1The Neutron Transport Equation is a linear, integro-diferential equationfor the neutron flux which depends of seven variables: three for position in space,two for directions, one for energy and one for time [2]. Solving this equation isa complex task, for which there are two possible approaches: deterministic andstochastic methods. Deterministic methods resort to the discretization of the trans-port equation with respect to its variables and converting the problem into a systemof algebraic equations to be solved. One of the main disadvantages of these methodsis the accuracy of its results, due to the discretization of the phase space (mesh orgrid resolution). On the other hand, stochastic methods simulate the physical trans-port problem randomly sampling the physical interaction of neutrons in a materialaccording to its reaction cross sections. Observables such as neutron flux, reactionrates, currents, among others are obtained by the expected value of N realizationsof the random sampling. The advantage of this method lies in the fact that does notresort to any approximation or discretization; its disadvantage is that the associatedstatistical uncertainty converges slowly as 1 / √ N , with N the number of particlessimulated. In this thesis, stochastic Monte Carlo method was used to solve the Neu-tron Transport Equation in fissile systems, approaching to be used in a completenuclear reactor model.While Monte Carlo methods are widely used in criticality and fixed sourcecalculations, where the system is supposed to be in stationary state, only recentlythere have been studies to include time dependence in neutron transport, takingadvantage of the better computing capabilities available. Some examples of thesestudies are the work of Snejitzer [3], Mylonakis [4] and Faucher [5], all of themfocused in the inclusion of time dependence together with the coupling of feedbackfrom thermal-hydraulics calculations.These investigations have in common the use of the customary group struc-ture for all the precursors. Each precursor group, which contains a number of differ-ent isotopes, is characterized by a grouped: i) decay constant; ii) relative yield; andiii) energy spectrum for the delayed neutron emission. This structure was proposedin 1957 by Keepin [6], and it is based in the assumption that the decay of the delayedneutron activity can be represented by a linear superposition of exponential decayperiods. Although this grouping is routinely used when performing deterministic orMonte Carlo simulations, it limits the possibility of studying the effect of changesin quantities such as the time evolution of the neutron flux, stimulated by new andimproved nuclear data on individual precursors.There has been a renewed interest in the measurement of nuclear decayproperties of the most neutron-rich nuclei, such as decay half-lifes, neutron emissionprobabilities and production yields [7], along with efforts from the InternationalAtomic Energy Agency Coordinated Research Project on a Reference Database for β -delayed Neutron Emission [8]. This scenario brings the opportunity to explorehow the new individual precursor data impacts on simulations of fissile systemsindividually or in different precursor groupings.The objective of this work is to explicitly include the time dependencerelated to the β -delayed neutron emission from individual precursors in a MonteCarlo simulation. This entails two challenges: to simulate the delayed emissionfrom precursors, and the inclusion of individual precursor data in the simulation.To include these modifications, the open source Monte Carlo code OpenMC waschosen [9].This work is divided in 5 chapters and 5 appendices. In Chapter 2, thetheoretical framework behind this work is presented. In particular, the NeutronTransport Equation (NTE) is examined, including the k -eigenvalue form and thepoint kinetics equation approximation. After that, the main features and differencesbetween prompt and β -delayed neutrons are discussed, along with the importantrole of β -delayed emission for nuclear reactor operation. The N -group structure fordelayed neutron precursors is also examined. Afterwards, the nuclear parametersneeded to understand the β -delayed neutron emission from individual precursors aredescribed, together with the nuclear data libraries used in this work, JEFF-3.1.1 [10]and ENDF/B-VIII.0 [1]. Finally, two of the approaches used to solve the NTEare discussed: deterministic and Monte Carlo methods. Related to the latter, adescription of variance reduction techniques is presented.In Chapter 3, methods used and developed to include the β -delayed neutronemission from individual precursors in a Monte Carlo simulation are discussed. Thefirst point addressed is about OpenMC, the code chosen to include the modificationsneeded to achieve the objectives of this work. This modified code will be known asTime-dependent OpenMC or OpenMC(TD). Afterwards, details on the methodologyto include time dependence in a Monte Carlo simulation are explained. The nextpart shows that time delay of β -delayed neutron emission, in an analog Monte Carlosimulation, entails large variance in the results obtained, due to the different timescales between the emission of prompt and delayed neutrons. To solve this prob-lem, forced decay of precursors is implemented in OpenMC(TD), but this strategyrequires population control of the neutron and precursor population. Regarding theinclusion of individual precursors, the steps taken to include them in OpenMC(TD)are: defining a precursor importance, so in the event of delayed neutron emission inthe simulation, it can be chosen which precursor will decay. This decay will have itsrespective precursor decay constant associated and the corresponding delayed neu-tron energy will be the average energy from the precursor delayed neutron spectrum.In Chapter 4, first the OpenMC(TD) code is tested in the context of time de-pendence and inclusion of individual precursors. With the tests successfully passed,OpenMC(TD) is used to obtain the neutron flux as a function of time in differentsystems, with different configurations and using different precursor structures. Thefirst system studied was a monoenergetic fissile system with 1-group precursor struc-ture, in subcritical, critical and reactivity insertion configurations. Afterwards, anenergy dependent, unmoderated U system was studied. This case was no longermonoenergetic, but energy dependent, using cross sections from JEFF-3.1.1 nucleardatabase. Two configurations were considered, subcritical and supercritical, and foreach the β -delayed neutron energies simulated were JEFF-3.1.1 and ENDF-B/VIII.0databases, in 1-group, 6-group, 8-group and 50 individual precursor structures. Thelast part of this chapter was related to simulations conducted in a light-water mod-erated, energy dependent system in a critical configuration with β -delayed neutronemission from 6-group, 50 individual, and 40 individual precursors.Finally, in Chapter 5 the conclusions and future perspectives of this workare presented. hapter 2Theoretical Framework In this chapter the theoretical framework behind this work is summarized. InSection 2.1 the Neutron Transport Equation (NTE) is examined, including the k -eigenvalue form (see Sec. 2.1.3) and the point kinetic equation approximation (seeSec. 2.1.4). Then, in Section 2.2 characteristics and differences between prompt and β -delayed neutrons are described, along with their important role for nuclear reactoroperation and the N -group structure for delayed neutron precursors (see Sec. 2.2.1).Afterwards, in Section 2.3 the nuclear parameters needed to describe the β -delayedneutron emission from individual precursors are shown, including the nuclear datalibraries used in this work (See Sec. 2.3.2). Finally, in Section 2.4 two of the ap-proaches used to solve the NTE are discussed, namely, deterministic (See Sec. 2.4.1)and Monte Carlo methods (See Sec. 2.4.2), where also a description of variance re-duction techniques is presented (See Sec. 2.4.3). Emphasis is given to time dependentphenomena, which are central to the challenges met throughout this work. The determination of the neutron distribution is the main problem of nuclear reactortheory because it determines the rate at which several nuclear reactions occur within5a fissile system. The knowledge about the neutron distribution also gives informationabout the stability of the fission chain reaction. The most general equation thatgoverns the process of neutron transport through a medium, this is, the motion ofneutrons as they stream through a system, is the Neutron Transport Equation (NTE)equation [2] (cid:20) v ∂∂t + ˆΩ · ∇ + Σ tot ( r , E ) (cid:21) ψ ( r , E, ˆΩ , t )= (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) Σ S ( E (cid:48) → E, ˆΩ (cid:48) → ˆΩ ) ψ ( r , E (cid:48) , ˆΩ (cid:48) , t ) + S ( r , E, ˆΩ , t ) . (2.1)In this equation the quantity to be determined is the neutron flux ψ , with E theneutron kinetic energy, ˆΩ the flux angular direction and v is an average neutronspeed. In this equation Σ tot and Σ s are the macroscopic total and scattering crosssections, respectively. Any external neutron source, such as fission neutrons, arerepresented by S . When including fission neutrons explicitly in Eq. (2.1), two contributions to thesource should be accounted for. The first term accounts for the prompt neutronsproduced in the nuclear fission process and is given by χ p ( E ) (cid:88) i (1 − β i ) (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) ν ( E (cid:48) )Σ f ( r , E (cid:48) ) ψ ( r , E (cid:48) , ˆΩ (cid:48) , t ) . (2.2)Here, Σ f ( r , E ) is the macroscopic fission cross section, ν ( E ) is the average number ofneutrons produced per fission, β i is the effective delayed neutron fraction per precur-sor group i , and χ p ( E ) the fast fission neutron spectrum. The second contribution tothe source term accounts for the delayed neutrons produced after the fission reactionand it reads (cid:88) l χ l ( E ) λ l C l ( r , t ) , (2.3)where the precursors are grouped in l groups according to their decay constant λ l , C l ( r , t ) represents the l -th precursor concentration and χ l ( E ) is the delayed neutronenergy spectrum for the l group. The precursor concentration, C l ( r , t ), changes intime as ∂∂t C l ( r , t ) = (cid:88) i β il (cid:90) dE (cid:48) (cid:90) d ˆΩ ν ( E (cid:48) )Σ f ( r , E (cid:48) ) ψ ( r , E (cid:48) , ˆΩ (cid:48) , t ) − λ l C l ( r , t ) , (2.4)where the first term on the right hand of Equation (2.4) stands for the producedprecursors while the left hand of the equation stands for decayed precursors. TakingEqs. (2.2), (2.3) and (2.4) into account, Eq. (2.1) for the neutron flux is reducedto [11] (cid:20) v ∂∂t + ˆΩ · ∇ + Σ tot ( r , E ) (cid:21) ψ ( r , E, ˆΩ , t )= (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) Σ S ( E (cid:48) → E, ˆΩ (cid:48) → ˆΩ ) ψ ( r , E (cid:48) , ˆΩ (cid:48) , t )+ χ p ( E ) (cid:88) i (1 − β i ) (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) ν ( E (cid:48) )Σ f ( r , E (cid:48) ) ψ ( r , E (cid:48) , ˆΩ (cid:48) , t )+ (cid:88) l χ l ( E ) λ l C l ( r , t ) . (2.5)It is important to remark that this equation relies on some assumptions: (i) neutronsare point-like; (ii) between two collisions neutrons travel in a straight line; (iii)neutrons do not interact with each other; (iv) collisions are instantaneous; and (v)materials do not change in time. Since the NTE features derivatives, appropriateinitial and boundary conditions must be specified for the neutron flux. The initialcondition can be the specification of the initial value for the neutron flux for allpositions, energies and directions: ψ ( r , E, ˆΩ ,
0) = ψ ( r , E, ˆΩ ) (2.6)The boundary condition will depend on the problem being studied, but usually theboundary conditions are: i) vacuum boundary condition; ii) reflective boundarycondition; and iii) known surface source. One of the most useful notations of the NTE is the steady-state form associatedwith the criticality of the system. In this problem the objective is the determinationof the k -eigenvalue ( k eff ) that shows the state of neutron multiplication in a fissilesystem or nuclear reactor. If k eff = 1, then the system is said to be in a criticalstate, if k eff < k eff > (cid:104) ˆΩ · ∇ + Σ tot ( r , E ) (cid:105) ψ ( r , E, ˆΩ )= (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) Σ S ( E (cid:48) → E, ˆΩ (cid:48) → ˆΩ ) ψ ( r , E (cid:48) , ˆΩ (cid:48) )+ χ p ( E ) (cid:88) i (1 − β i ) (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) ν ( E (cid:48) )Σ f ( r , E (cid:48) ) ψ ( r , E (cid:48) , ˆΩ (cid:48) )+ (cid:88) l χ l ( E ) (cid:88) i β il (cid:90) dE (cid:48) (cid:90) d ˆΩ ν ( E (cid:48) )Σ f ( r , E (cid:48) ) ψ ( r , E (cid:48) , ˆΩ (cid:48) ) . (2.7)By defining the net disappearance operator L as Lf = ˆΩ ·∇ f +Σ tot ( r , E ) f − (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) Σ S ( E (cid:48) → E, ˆΩ (cid:48) → ˆΩ ) f ( r , E (cid:48) , ˆΩ (cid:48) ) , (2.8)and the total fission operator F as F f = χ p ( E ) (cid:88) i (1 − β i ) (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) ν ( E (cid:48) )Σ f ( r , E (cid:48) ) f ( r , E (cid:48) , ˆΩ (cid:48) )+ (cid:88) l χ l ( E ) (cid:88) i β il (cid:90) dE (cid:48) (cid:90) d ˆΩ ν ( E (cid:48) )Σ f ( r , E (cid:48) ) f ( r , E (cid:48) , ˆΩ (cid:48) ) , (2.9)Eq. (2.7) can be rewritten as Lψ ( r , E (cid:48) , ˆΩ (cid:48) ) = F ψ ( r , E (cid:48) , ˆΩ (cid:48) ) . (2.10)By imposing that the system should be critical, the k -eigenmodes can be found L ψ k ( r , E (cid:48) , ˆΩ (cid:48) ) = 1 k eff F ψ k ( r , E (cid:48) , ˆΩ (cid:48) ) . (2.11) One theoretical approximation used to study the transient behaviour of a nuclearreactor is the point kinetic approximation [12]. In this case, the flux is assumed to bea separable function of space and time and the equations are obtained by weightingthe transport equation by the adjoint flux.To obtain the pertinent equations first Eq. (2.5) can be written as1 v ∂∂t ψ ( r , E, ˆΩ , t ) + L ψ ( r , E, ˆΩ , t ) = F p ψ ( r , E, ˆΩ , t ) + (cid:88) l χ l ( E ) λ l C l ( r , t ) , (2.12)where the prompt fission operator, F p , is defined as F p f = χ p ( E ) (cid:88) i (1 − β i ) (cid:90) ∞ dE (cid:48) (cid:90) π d ˆΩ (cid:48) ν ( E (cid:48) )Σ f ( r , E (cid:48) ) ψ ( r , E (cid:48) , ˆΩ (cid:48) , t ) f. (2.13)The adjoint equation to the k -eigenmodes Equation (2.11) is L † ψ † k ( r , E, ˆΩ ) = 1 k eff F † ψ † k ( r , E, ˆΩ ) , (2.14)where ψ † k is the adjoint eigenmode for the neutron flux and the L † is the adjoint ofthe operator.To derive the point kinetic equations, the transport Eq. (2.12) is multipliedby ψ † k and Eq. (2.14) is multiplied by the neutron flux ψ ( r , E, ˆΩ , t ). The resultingequations are integrated over space, energy and angle and then are substracted fromeach other, obtaining ∂∂t (cid:104) ψ † k , v ψ (cid:105) = k − k (cid:104) ψ † k , F ψ (cid:105) − (cid:104) ψ † k , F d ψ (cid:105) + (cid:88) l λ l (cid:104) ψ † k , χ ld C l (cid:105) . (2.15)It is assumed that the neutron flux can be factorized as an amplitude factor thatonly depends on time and a time-independent flux shape factor: ψ ( r , E, ˆΩ , t ) = s ( t ) ψ k ( r , E, ˆΩ ) , (2.16)where ψ k ( r , E, ˆΩ ) is the fundamental k eigenmode and s ( t ) is an amplitude factorthat only depends on time. Thus, Eq. (2.15) the amplitude of the neutron populationsatisfies ∂∂t s ( t ) = ρ − β eff Λ eff s ( t ) + (cid:88) l λ l c l ( t ) , (2.17)where ρ is the reactivity given by ρ = k − k , (2.18)0the effective delayed neutron fraction is β eff = (cid:104) ψ † k , F d ψ k (cid:105)(cid:104) ψ † k , F ψ k (cid:105) , (2.19)the effective mean generation time isΛ eff = (cid:104) ψ † k , v ψ k (cid:105)(cid:104) ψ † k , F ψ k (cid:105) , (2.20)and the effective precursor concentration c l ( t ) = (cid:104) ψ † k , χ ld C l (cid:105)(cid:104) ψ † k , v ψ k (cid:105) . (2.21)By proceeding in a similar way, an equation for the precursor concentration can bederived, which is coupled to Eq. (2.17), ∂∂t c l ( t ) = β l Λ eff s ( t ) − λ l c l ( t ) , (2.22)where β l is the effective delayed neutron fraction for the precursor family l .Parameters β eff and Λ eff are called effective because they have been weightedby the adjoint flux ψ † k , which can be interpreted as the neutron importance. Physi-cally, the neutron importance at a given point in phase space is proportional to theasymptotic neutron population of an hypotetical neutron introduced into a criticalreactor at the same point in phase space. In section 2.1 the transport equation was presented. For nuclear fission present ina fissile system, the source term is comprised by two terms, the prompt and the delayed fission term . In fission events, two types of neutrons are released: promptand delayed neutrons.Prompt neutrons are released almost instantaneously ( ∼ − s) after fis-sion and are emitted with an average energy of 2 MeV [13]. Since the fission cross1 - - - - -
10 1 10 Energy (MeV) C r o ss - s e c t i on ( b ) uranium-235 Figure 2.1: Fission Cross sections for
U. Data was retrieved from the JEFF-3 . . ν p .On the other hand, delayed neutrons are emitted between 10 − s and 10 safter a nuclear fission event and made up about 1 % of the total neutrons releasedduring fission. This delayed emission is produced when a fission product ( Z, N )decays through a β − process with a father-daughter mass difference Q β . If Q β isgreater than the neutron separation energy S n , then excited states in the daughternucleus ( Z + 1 , N −
1) can be populated. This nucleus can in turn decay to thenucleus ( Z + 1 , N − β − decay corresponding to the( Z, N ) nucleus. This parent β -decay nucleus ( Z, N ) is known as delayed neutronprecursor or precursor . To illustrate this process, in Fig. 2.2 the decay scheme of theprecursor Br is shown. In this scheme it can be seen that Br can decay through β − to a state in Kr ∗ , followed by the subsequent decay of Kr ∗ to a state in Krvia neutron emission. The delayed time of this process is given by the parent halflife, which is 55 . ν d . The fraction of total fission delayed neutrons is denoted by β and isdefined as β = ¯ ν d ¯ ν . (2.23)2Figure 2.2: Decay of the Br delayed neutron precursor .In a fission chain reaction a large number (about 270 for
U) of delayedneutron precursor isotopes can be produced [14]. It has been customary to groupthese precursors in 6 groups, characterized by its half-lives. Each precursor groupcontains a number of different isotopes. These groups are described in Table 2.1for
U, where for each group is shown: i) the group mean energy, ¯ E , ii) its half-life, T / and, iii) its relative yield, given by β i /β , where β i is the delayed fractionconsidering only the i -th group and the total delayed delayed fraction is (cid:80) i β i = β .This group structure was proposed by Keepin [6], who assumed that the decay ofdelayed neutron activity with time could be represented by a linear superposition ofexponential decay periods. Group Precursors ¯ E (MeV) T / (s) β i /β Br,
Cs 0 .
40 54 .
51 0 . I, Br 0 .
47 21 .
84 0 . I, Br, Rb, Rb 0 .
44 6 .
00 0 . I, Xe, Kr, Kr, Br, Br 0 .
55 2 .
23 0 . I, Cs 0 .
51 0 .
496 0 . .
54 0 .
179 0 . Table 2.1: 6 delayed neutron precursor groups for
U fission. In 6-th group severalisotopes of Br, As, Rb are included.3
Neutron Fissilematerial Fission Fission productPrecursorPrompt neutrons EmitterDelayed neutron Stable nucleus
Figure 2.3: Schematic representation of the prompt and β -delayed neutron emission. Delayed neutrons are important in the operation of a nuclear reactor due to thetime delay that they introduce to the system, needed to control the state of thereactor through mechanical means, such as control rods. To illustrate this point, thevariation of the neutron density without taking into account delayed neutrons is dndt = k − (cid:96) n ( t ) = ∆ k(cid:96) n ( t ) , (2.24)where n ( t ) is the neutron density, k is the multiplication factor and (cid:96) is the promptneutron lifetime, taken as the average time a neutron stays in the system beforeleaking or being absorbed. The solution to this Eq. is n ( t ) = n exp (cid:18) ∆ k(cid:96) t (cid:19) ≡ n exp (cid:18) tτ (cid:19) , (2.25)with n the initial neutron density. The rate at which the reactor power increasesis given by ∆ k/(cid:96) . The reciprocal of this quantity is the reactor period, τ = (cid:96)/ ∆ k ,namely the time needed for the reactor power to grow by a factor of e .If there are no delayed neutrons, then the mean neutron lifetime is the meanprompt neutron lifetime (i.e. (cid:96) = (cid:96) p ) which in a light-water reactor is about 10 − s [2].4Then, if there is a positive reactivity insertion of 10 pcm , the multiplication factorwould change from k = 1 . k = 1 . k = 0 . .
01% of k ) andthe reactor period is τ = 10 − / . . ≈ − β eff ) of the totalneutrons have the prompt neutron lifetime (cid:96) p , while the delayed neutrons, a fractionof β eff the neutrons, live longer with a lifetime ( T avg + (cid:96) p ). This implies that theneutron lifetime in this case is (cid:96) = (cid:96) p (1 − β eff ) + ( T avg + (cid:96) p ) β eff ≈ β eff T avg , (2.26)and if Eq. (2.26) is used in Eq. (2.25) the reactor period in this case would be τ = 0 . / . ≈
100 s. With this reactor period, in a second the power would riseby a factor of ≈ .
1, making easier to control the behavior of the reactor in the faceof reactivity insertions.Then, due to the effect of delayed neutrons, the period of the reactor in-creases and the rise in reactor power slows down, making possible to control thereactor by mechanical means.
In this section the quantities of interest needed to characterize individual precursornuclides are presented. Then, the two nuclear databases considered in this work,ENDF/B-VIII.0 and JEFF-3 . .
1, are briefly described along with some of the differ-ences found in the course of this work.
Although the 6- or 8- group structure is widely used in reactor calculations [15],nowadays there is not a solid reason to keep this aggregation. More and better nuclear Per-cent mille: one-thousandth of a percent.
F Y ), the precursor decay constant( λ ) and the precursor delayed neutron emission probability ( P n ). To study theeffect of individual precursors in a Monte Carlo simulation, information related tothe energy distribution of the emitted neutrons is also needed. Now each of thisquantities will be reviewed [17]: Fission Yield (
F Y ) : Two yields can be defined, one is the Independent FissionYield (IY) , which is the average number of atoms of a specified nucleus produced byone fission, after the emission of prompt neutrons and excluding radioactive decay.For example, Br is one of the most prominent delayed neutron precursor and its IY is 0 . Br are produced. The other is the
Cumulative Fission Yield(CY) , which is the number of atoms of a specific nuclide produced directly and viadecay of precursors per one fission reaction.
Precursor Decay Constant ( λ ) : The decay constant represents the probability fora nucleus to decay per time unit. The decay probability of a precursor is proportionalto the number of nuclei. λ = − dN/dtN (2.27) Precursor delayed neutron emission probability ( P n ) : For β -delayed neutronemission to occur, the β − decay energy ( Q β ) must be larger than the neutron sep-aration energy ( S n ) of the decay daughter. The precursor delayed neutron emissionprobability represents the probability of one or more neutron emission. Precursor delayed neutron spectrum : the energy distribution of neutrons emit-ted by each precursor is characterized by its spectrum. At this moment, the ENDF/B-VIII.0 database has only 34 evaluated experimental spectra while the others comefrom QRPA calculations [14]. In this work the “mean energy” of the delayed neutronsemitted was used.6
Average delayed neutron yield ( ν d ) : Also known as the average number ofdelayed neutrons produced per fission, it can be either measured [6], or calculatedusing the cumulative yield and the precursor delayed neutron emission probability, ν d = N (cid:88) i CY i P n,i (2.28)where, N is the number of precursors. In 1990 the Nuclear Energy Agency establisheda Working Party on International Nuclear Data Evaluation Co-operation (WPEC)to “promote the exchange of information on nuclear data evaluations, validation andrelated topics. Its aim is also to provide a framework for co-operative activitiesbetween members of the major nuclear data evaluation projects”. The Subgroup 6(WPEC-6) in particular had the objective of reducing the uncertainties on delayednuclear data and in this context, the recommended value for the average delayedneutron yield for U thermal is 1 . × − [15]. A nuclear data library is a dataset of stored nuclear data in a certain format [18].Nuclear data derived from the combination of experimental data and nuclear physicsmodels are known as evaluated nuclear data libraries. The standard format for thestorage of nuclear data is the ENDF-6 format (
Evaluated Nuclear Data File ) [19]. Inthe course of this work nuclear data from two libraries was studied, the
JEFF and
ENDF/B libraries. The JEFF (
Joint Evaluated Fission and Fusion File ) library iscreated by the OECD/NEA and the version used in this thesis was JEFF-3 . . Evaluated Nuclear Data File / B ) library is created by the
Cross Sec-tion Evaluation Working Group and the version used was ENDF/B-VIII.0 [1]. Bothlibraries contain radioactive decay data sub-libraries, where the yields, branchingratios and delayed neutron spectra can be found.It was found that the data from both libraries do not agree with each otherand that there are important differences between quantities of interest. To show this,the independent yield, cumulative yield and branching ratio is shown for some of themain precursors from the 6-group structure [20], as it is shown in Table 2.2.A summary of the differences between both libraries is shown in Table 2.3.7 IY (% per fission) CY (% per fission) β n (% per fission) Group Nucleus
ENDF/B-VIII.0 JEFF-3 . . . Br 1 . . . . . . Br 1 . . . . . . I 2 . . . . . . Br 1 . . . . . . Rb 3 . . . . . I 1 . . . . . . As 0 . . . . . Br 0 . . . . . . Rb 1 . . . . . . I 0 . . . . . . Br 0 . . . . Rb 0 . . . . . . Rb 0 . . . . . . Table 2.2: Independent yields, cumulative yields and branching ratio values foundin JEFF 3 . Br the difference inthe values for the independent yield is the order of 10% and in extreme cases thisdifference can take values up to 170%, as in the case of As.
Group Nucleus ∆ IY (%) ∆ CY (%) ∆ β n (%)1 Br 10 5 − . Br 6 2 1 . I 11 14 − . Br 20 20 2 . Rb 8 1 − I 3 1 − . As 14 35 − Br 14 14 − . Rb 11 9 − I 24 23 − Br 32 32 100 Rb 15 15 − . Rb 60 50 0 . Table 2.3: Independent yields, cumulative yields and branching ratio values foundin JEFF 3 . . CY ’sused were taken from JEFF 3 . .
1, while the P n were taken from ENDF-B/VIII.0.This pairing is the recommended when comparing the ν d calculated using the sum-mation method given by Eq. (2.28) with the experimental value [17].8 Basically there exists two different approaches to solve the Neutron Transport Equa-tion: deterministic and stochastic (Monte Carlo) techniques. Although they aim tostudy the same physical problem, they are different in their approach and techniquesused to solve the problem. They have complementary advantages and disadvantages.
Deterministic methods model the physical problem ignoring the random aspect ofindividual particle histories and then they solve the NTE by discretizing these equa-tions with respect to each of its variables and converting the problem into a systemof algebraic equations that has to be solved [11]. One of the strategies used to cal-culate the neutron flux use the quasi-static method, developed in the 1950s [21].This method resorts to an approximation where the flux is factored as a productbetween a shape function and an amplitude function. The shape function can beobtained through stationary state calculations, using discrete ordinates [22, 23] orMonte Carlo [24,25]. To obtain the time-dependent amplitude function, diffusion [26]or S n methods [27] can be used. One feature of all of these methods is that theydiscretize the phase space: difussion theory assumes that neutrons diffuses throughthe medium following Fick’s law and ignores the angular dependence of the flux.More advanced methods such as the S n method does take into account the angulardependence of the flux, but this dependence is discretized and neutrons are trans-ported though discrete angles. With the use of these techniques it can be possibleto refine the modelling of the angular dependence of the flux, but it can be complexto use the necessary number of angles to obtain a good solution for the flux. One ofthe main disadvantages of these methods is the constraints in the resolution of thediscretization grid, since memory is required to store the unknown variables. Witha coarser grid, higher discretization errors are obtained, so limitations in memorylimit the accuracy of deterministic methods [28].9 Unlike the methods described in the preceding section, the Monte Carlo method(which is an stochastic method) do not solve the Neutron Transport Equation ex-plicitly, but simulate the physical problem by transporting the neutrons through themedium. The physical processes involved in the evolution of the neutron popula-tion is governed by probability distributions. In the application of the Monte Carlomethod to neutron transport, a stochastic model is simulated, and then the expectedvalue of some random variable is equivalent to the value of a physical quantity thatis to be determined. This quantity is estimated using the average of independentsamples that represent the random variable.To illustrate this point, the procedure to carry out a Monte Carlo simulationwill be outlined [28]. For simplicity a time-independent fixed source problem in ahomogeneous medium will be considered. In this problem a source and a detector inthe phase space must be simulated, and the detector response will be the quantityto be estimated, this is, the contributions of neutrons reaching the detector will becollected. The idea is to simulate N neutrons, sampling the source distribution tofind initial energy, position and direction for each neutron. The emitted neutronsare then transported. The distance d that each neutron of energy E travels betweentwo interactions is exponentially distributed and given by d = − ln( ξ )Σ t ( E ) , (2.29)where ξ is a random number, sampled from a uniform distribution between [0 , t is the total macroscopic cross section. If d is larger than the distance to theboundary of the next volume, then the particle is stopped at that boundary and anew path is sampled using Eq. (2.29). At the new position the interacting nucleus i needs to be sampled, which will be chosen with probability p i = Σ t,i ( E )Σ t ( E ) , (2.30)where Σ t,i ( E ) is the total cross section for nucleus i . Once the interacting nucleus issampled, the specific interaction occurs with a probability p i,x = σ i,x σ i,t , (2.31)0where σ i,x is the microscopic cross section for the interaction x and nucleus i .After the interaction the neutron can be eliminated if absorbed or if it leavesthe simulation world. Otherwise a new path is sampled and the process starts again.Neutron contributions are accumulated when they reach the detector. After the N particles are transported the process is repeated M times with a different randomseed each time.The tallies collected are averaged over the M experiments. The associateduncertainty is calculated using the variance and its inversely proportional to thesquare root of the number of particles simulated. This means that the uncertainty ofthe result obtained in a simulation can be improved by simulating a larger numberof particles. Although the number of particles required is problem dependent, itis usually quite large and this implies that Monte Carlo simulations are very timeconsuming. Fortunately, Monte Carlo algorithms are specially suited for parallelcomputing [29], which allows to speed up, in principle, by the order of the processoravailables. The idea is that each processor simulates its own number of particles, andwhen each processor have completed the transport, the final results are collected. A Monte Carlo simulation as described in section 2.4.2 requires knowledge of theprobability distribution that governs the physical process that is used to calculatethe expected value. In other words, in this method the computation describes how aparticle would behave in an equivalent physical experiment. This method is knownas analog
Monte Carlo simulation [11]. There are some experimental setups where,for example, the detector counting rate could be too low or, for a shielding problem,there are too few initial particles that reach the region of interest. In those cases,longer detection times or several repetitions of the experiment might be necessaryto achieve an acceptable uncertainty. If one of these physical systems would besimulated using Monte Carlo, large number of particles would be required in orderto achieve a reliable estimate of the quantity being studied. But the simulation timeis governed by the number of particles simulated, which means that the simulationwould require very long computation times. One way to overcome this problem is1through a non-analog
Monte Carlo simulation [30], which is a modification of the analog
Monte Carlo simulation where the physical probability distribution is modifiedin order to promote the occurrence of a given event (for example, to make that moreparticles can reach the detector). To keep the results unbiased, a compensationhas to be applied elsewhere. For this purpose a statistical weight is defined andassigned to each particle at the beginning of the simulation. Then, this weightcan evolve along the simulation to counterbalance the changes that occur when theprobability of a physical process is altered. One of the variance reduction techniquesused in this work is survival biasing , which must be used in conjunction with apopulation control technique called
Russian roulette . In survival biasing (also knownas implicit absorption ), absorption reactions are prohibited to occur and instead atevery collision the statistical weight of the particle, w new is reduced by the probabilitythat the absorption occurs: w new = w (cid:18) − σ a ( E ) σ t ( E ) (cid:19) , (2.32)where σ a ( E ) and σ t ( E ) are the absorption and total microscopic cross sections, re-spectively, and w is the statistical weight of the particle before the collision. It isimportant to notice that survival biasing can reduce the weights of the particles tovery low values. In that case, particles of low statistical value slow down the calcu-lation, while contributing very little to the statistics. This means that this methodmust be combined with another method capable of stochastically killing particles.This method is called Russian roulette . If a particle falls below some thresholdweight, then a random number is generated. If the random number is below theinitial weight, then the particle is killed . Otherwise, it survives and its weight is setto some predefined value. In the context of Monte Carlo simulations, to kill a particle is to remove it from the simulation. hapter 3Methodology
In this chapter, the methodology used to include the β -delayed neutron emissionfrom individual precursors in transient Monte Carlo simulation is discussed. Firstly,in Section 3.1 the Monte Carlo OpenMC code is described, along with explanationwhy it becomes suitable for this work, and a benchmark calculation result is pre-sented (See Sec. 3.1.1). After that, in Section 3.2 a discussion on how the timedependence is treated. Following, Section 3.3 addresses precursors, including conse-quences of β -delayed neutron emission in the context of a Monte Carlo simulation(see Sec. 3.3.1). This comprises how individual precursors are implemented in thecode (see Sec. 3.3.2), and the strategy to overcome large variances associated withthe different time scales between prompt and delayed neutrons (see Sec. 3.3.4). Af-terwards, in Section 3.4 the issue of how to sample a proper initial source to starta Monte Carlo transient simulation is discussed. Finally, in Section 3.5 the methodchosen to enforce population control is described. The OpenMC [9] code is relatively new, an open-source code for particle transportdeveloped at the Massachusetts Institute of Technology in 2013. This code is capableof simulating neutrons in fixed source, k -eigenvalue, and subcritical multiplicationproblems. The geometry is built using a constructive solid geometry. The code223supports both continuous-energy and multigroup transport. The continuous-energynuclear cross section data follows the HDF5 format [31] and is generated from ACEfiles produced by NJOY [32]. Since this code is open source, its use is not subject tolicensing, with no restrictions on modifications, developments and addition of newcapabilities. As a first step before beginning the development of new capabilities for OpenMC,benchmark calculations were performed to further validate the code. In order to dothis, and during the author’s first doctoral internship at the Bariloche Atomic Cen-ter, the OpenMC code was used to model and calculate the Effective MultiplicationFactor of the RA-6 research reactor. The result obtained was compared with the ex-perimental values from the ICSBEP
International Handbook of Evaluated CriticalitySafety Benchmark Experiments [33, 34], and with the result obtained when modelingthe reactor using the Monte Carlo transport code MCNP [35].Another parameter that can be calculated is the Effective Delayed NeutronFraction, β eff , which was mentioned in Sec. 2.1.4. Formally, the adjoint neutron fluxis required to calculate this parameter, but it can be estimated using the promptmethod [36]. This method assumes that the value of β eff is given by β eff ∼ − k p k eff , (3.1)where k p is the effective multiplication factor obtained from a criticality calculation,but without taking into account the contribution from β -delayed neutron emission.The advantage of this method is that the adjoint flux is not needed to calculate theeffective delayed neutron fraction. Thus, the capability to run a criticality calculationwithout delayed neutrons was added to OpenMC, which enabled the estimation of β eff in two steps. In MCNP6 the same feature can be achieved by using the TOTNU NO card to perform a criticality calculation only with prompt neutrons.
Description of the RA- reactor The RA-6 (Spanish acronym for
Argentina Reactor, Number 6 ) is an open4pool research reactor with a nominal power of 3 MW, located at Bariloche AtomicCenter, a nuclear research center in San Carlos de Bariloche, R´ıo Negro, Argentina.The core of the reactor is made up of an array of flat plates MTR-type fuel elementswith 20 % enriched uranium located inside a stainless steel tank filled with deminer-alized water that acts as a coolant, moderator, reflector and shielding in the axialdirection. Four Ag-In-Cd absorber elements are the control elements. The modelwas the one included for the ICSBEP benchmark evaluation, with added graphitereflectors [34] and –since in Monte Carlo codes it is possible to model the reactorgeometry in detail– fuel elements were modeled explicitly, such as cadmium wires,water gaps, guides and nozzles. The model also included the supporting grid for thecore and BNCT filter.
Simulation parameters and results for k eff calculation In OpenMC and MCNP the criticality calculation was peformed using 8050batches , 50 skipped and 10000 particles per batch. The neutron cross sectiondatabase used was ENDF/B-VII.1. Results obtained are summarized in Table 3.1Magnitude OpenMC MCNP Benchmark k eff . . . k p . . − β eff As stated in the Introduction, the main objective of this thesis is to study the in-clusion of time-dependence in a Monte Carlo simulation, considering the delayedemission from the neutron precursors present in a fissile system. To this end severalissued must be addressed, which will be discussed in the remainder of this chapter. the total number of source particle simulated is broken up into a number of batches. skipped cycles will be discarded before data accumulation begins. In a stationary Monte Carlo transport simulation time is not explicitly present. Thefirst step to perform Monte Carlo kinetic simulations is to add a new label t to theparticles, serving as a clock with value updated using the kinetic energy and thedistance traveled by the neutron between events. This time is set to zero ( t = 0) atthe beginning of the simulation and is updated as the particle is transported in thesimulation. To simulate transient events in fissile systems the evolution was divided in discretetime intervals. There are two reasons for this: First, the variance reduction andpopulation control techniques require a time grid to be applied. The second reasonis that changes in the geometry or reactivity of the system can take place in atransient simulation, changes that can be introduced at the end of a time interval. Itis important to notice that the size of the time intervals can be choosen freely and theydo not affect the validity or accuracy of the results obtained from the simulation.When a particle crosses a time boundary, its trajectory is stopped exactly at theboundary, with the spacial position that corresponds to the time boundary, then theparticle is stored to continue the simulation at the next time step.
In order to tally the measured quantities in time, the tallies in OpenMC were modifiedand a new filter was added. This time filter added the capability to monitor the timeevolution of any of the tallies already present in the code.6
In this section the time delay of the β -delayed neutron emission and its consequencesin the context of a Monte Carlo simulation are explored. A key point of this discussionis the large variance in the simulation results if an analog Monte Carlo method wasused, issue that will be dealt with at the end of the section. Then the inclusion of β -delayed neutron emission from individual precursors, one of the main objectivesof this work, is addressed. Following this discussion, the precursor particle definedin this work, for the simulation was presented. In the last part of this section,techniques chosen to solve the problem caused by time differences between promptand delayed emission of neutrons are described. The delayed neutron precursor decay is a stochastic process, which can be describedby p i ( t ) = λ i e − λ i ( t − t ) θ ( t − t ) , (3.2)where p i ( t ) is the probability the i -th precursor decay at a given time t , λ i is the decayconstant, t is the time when the precursor was created and θ the Heavyside function.Given this probability, an analog Monte Carlo simulation could be performed to,in principle, describe what happens in a fissile system: at time t f of the fissionevent, ν p prompt neutrons are produced and then, at a time t = t f + t d , ν d delayedneutrons are inserted into the simulation. Time t d is sampled from Eq. (3.2) andthe energy is determined from the precursor delayed neutron energy distribution (seeSection 2.3.1).Although this strategy emulates what happens in a nuclear reactor, a largevariance in the results obtained due to the difference in the time scales associatedwith prompt and delayed events. Indeed, as it was mentioned in Section 2.2.1 thereexists a time delay between the nuclear fission event and the emission of a delayedneutron from the decay of a precursor. The average lifetime of a prompt neutronin a light water reactor is ∼ − s and the average length of a fission chain in asystem close to critical is ∼
150 neutrons [3]. This implies that the average lifetime7of a neutron chain is ∼ − s. At the same time, a prompt fission chain willproduce on average one precursor, which in turn will decay to a delayed neutronand then produce a new fission chain in a few seconds. During this time therewould be no new neutrons produced in an analog Monte Carlo simulation, as it isshown in Fig. 3.1. This lack of particles would in turn lead to large variance in thequantities scored. In an actual fissile system this does not happen because of thelarge number of neutrons produced so the effect is averaged out. Of course, due tolimitations imposed by computer calculation power and memory, it is not possibleto simulate this many fission chains. Due to this fact, and in order to obtain resultswith acceptable statistics, delay of precursors decay must be simulated in anotherway. Precursor
Delayedneutron N e u t r o n p o p u l a t i o n Precursor
Delayedneutron
Figure 3.1: Schematic representation of the time scales associated to the delayedneutron emission and the lifetime of the prompt chains. This different time scalesproduce large variance in the quantities scored.8
In 1957, Keepin measured the periods, relative abundances and yields of delayedneutrons from fission. He had the idea of grouping the β − delayed neutron emittersinto groups according to their half-lives and assuming that the total emission ratecould be represented as a sum of exponential functions [6]. It is important to notethat the number of groups is arbitrary, considering that the total number of precur-sors produced in U fission is more than 270. Nonetheless, Keepin found that a sixgroup representation properly fitted the measured experimental activity.At the time of the writing of this thesis, there were no published MonteCarlo codes for neutron transport in fissile systems which include the delayed neutronemission from individual precursors, i.e. all of the existing codes use the groupstructure to take into account β -delayed neutron emissions and insert a delayedneutron directly into the simulation. In the case of this work a precursor is created,and then this precursor can decay, emitting a delayed neutron which is inserted intothe system. To further illustrate this point, the steps needed to take into account β -delayed neutron emissions will be outlined using, i) the group structure or, ii) theindividual precursors. i) β -delayed neutron emission with group structure If a delayed fission is sampled, the next step is to choose which precursorgroup will be sampled. Here, the relative yield, is utilized. If the j -th group ischosen, then the delayed time associated with the delayed emission will be sampledusing the group decay constant and Eq. (3.2). Finally, the delayed neutron energywill be chosen from the j -th delayed neutron energy group spectra [1].9 ii) β -delayed neutron emission with individual precursors On the other hand, if a delayed fission is sampled and the delayed neutronemission from individual precursors is being simulated, instead of directly inserting adelayed neutron, a precursor is produced. The next step is to choose which precursornuclide will decay. In order to do this, the precursor importance or relative yield ofthe individual precursor i -th, I i , is defined as [17]: I i ≡ CY i P n,i ν d , (3.3)with CY i the cumulative fission yield, P n,i the precursor delayed neutron emissionprobability, and ν d the average delayed neutron yield. Once the precursor has beenchosen, the delayed time associated with this emission is sampled using the precursordecay constant and Eq. (3.2). Finally, the delayed neutron energy will be the averageenergy from the corresponding precursor delayed neutron spectrum.Another point to consider is the number of precursors to include in thesimulation, for which the precursor importance is useful because it shows the fractionof the total delayed neutron yield that the precursor represents (i.e. how importantit is). As an example, when using the cumulative yields from JEFF-3 . . β -delayed neutron emission probabilities from ENDF/B-VIII.0, the averagedelayed neutron yield obtained is 1 . × − . Then from the values presented inTable 2.2, the importance for any given precursor can be calculated. For example,for I, the precursor importance obtained is 16 . I decay account for 16 .
26% of the total β -delayedneutron emission.So, although there are data for 269 precursors, in this work only 50 willbe included in the simulation. To justify this choice the precursors were orderedby importance using Eq. (3.3) and then the cumulative importance ( (cid:80) i I i ), wascalculated. It was found that the first 50 precursors account for 99 .
16% of thetotal delayed neutron yield, which means that the remaining 219 precursors havea combined importance of 0 . β -delayed neutron emission using the group structure and the delayedneutron emission when using individual precursors. Quantity N-group structure This workRelative abundance β i /β with 1 < i < CY i P n,i ) /ν d with 1 < i < Decay constants
Precursors in 6- or 8- groups 50 individual precursors
Energy spectra
Precursors in 6- or 8- groups 50 individual precursors
Table 3.2: Summary of the differences when including the β -delayed neutron emissionusing the precursor group structure or the individual precursorsFinally, it is worth mentioning that the choice of including 50 out of the 269precursors was made taking into account the calculation time and the cumulativeimportance of these 50 precursors, but should need arise, the code developed canhandle the whole set of precursors. The first step to include the precursor decay in the simulation involves adding theprecursors in the simulation. So a new particle type is defined in the code, the precursor particle . All precursors (or precursor groups) are combined into a singleprecursor particle [3]. The decay probability for this particle is given by p combined ( t ) = (cid:88) i Γ i λ i e − λ i ( t − t ) θ ( t − t ) , (3.4)with, t the time when the precursor was created, and Γ i a factor that depends onwhether precursor groups or individual precursors are being considered: Γ i = β i β , for precursor group I i , for individual precursor . (3.5)Here β i is the delayed fraction for the i -th precursor and (cid:80) i β i = β , with β thetotal delayed neutron fraction. I i is the precursor importance for the i -th precursor.Fractions Γ i must be defined differently in some cases, as will be shown in Section 3.4.In principle, the statistical weight of a delayed neutron emitted from the β -decay ofa precursor is given by w d ( t ) = w c (cid:88) i Γ i λ i e − λ i ( t − t ) θ ( t − t ) , (3.6)1with w c the main precursor weight. This weight is the number of physical precursorsthat this precursor particle represents at the time of its creation in the simulation. Itmust be noted that this weight does not change with time and it can only be alteredby means of variance reduction techniques, as it will be explained in Section 3.5.After the precursor decay is produced, the energy of the emitted delayed neutronmust be chosen from the corresponding precursor or precursor group. The probabilityof choosing the i -th group or precursor is a function of time given by P i ( t ) = Γ i λ i e − λ i ( t − t ) (cid:80) i Γ i λ i e − λ i ( t − t ) . (3.7)This means that this probability must be evaluated at the time of decay to selectthe correct group or precursor for the energy spectrum.Aside from the main precursor weight w c , there is another statistical weightwhich will be utilized during this work. This is the weight of the precursor at atime t and it represents the number of physical precursors that a precursor particlerepresents at a given time t and is given by w p ( t ) = w c (cid:88) i Γ i e − λ i ( t − t ) . (3.8)The last statistical weight that can be utilized is the expected delayed neutron weight.The precursor interacts with the system through delayed neutrons, so the weight ofthe delayed neutrons can be used for variance reduction. The problem is that thedecay time is not known a priori , so this weight is defined as [3] w d,av = 1∆ t (cid:90) t +∆ tt w d ( t ) dt, (3.9)where t is the start of the next interval. Using Eq. (3.6), the expected delayedneutron weight becomes w d,av = w c (cid:88) i Γ i ( e − λ i ( t − t ) − e − λ i ( t +∆ t − t ) ) . (3.10) As explained in Section 3.3.1, a direct simulation of delayed neutron precursor decaywould lead to significant variance in the system, so that another way to simulate2the precursors must be utilized. Since this variance is caused by the fact that thereare too few fission chains per unit of time caused by delayed neutron decays, onestrategy would be to modify the precursor decay probability, forcing the decay of allthe precursors in each interval and thus having more delayed neutrons present. Inthis technique, called “forced decay” [37], the sampling of the delayed neutrons isbiased and the Monte Carlo fair game is preserved by altering the statistical weight ofthe emitted delayed neutrons. Regarding the biased decay probability, the simplestchoice would be a uniform decay probability, forcing the decay of all of the precursorsin each one of the time intervals defined in Section 3.2.2. With this choice the biaseddecay probability is ¯ p ( t ) = 1 t j +1 − t j = 1∆ t , (3.11)where t is the time when the forced decay happens and ∆ t is the size of the timebin. To ensure an unbiased result the weight of the delayed neutrons produced byforced decay is adjusted to w d ( t ) = p ( t )¯ p ( t ) = w c ∆ t (cid:88) i Γ i λ i e − λ i ( t − t ) with t j < t < t j +1 (3.12)where w c is the statistical weight of the precursor. The delayed neutron producedwill be transported and may in turn cause new fissions. Once the delayed neutronof weight w d ( t ) has been created during the corresponding time interval between t j and t j +1 , the precursor is not eliminated from the simulation. Instead, it is added toa precursor bank with weight w p ( t ) = w c (cid:88) i Γ i e − λ i ( t i +1 − t i ) , (3.13)where it will undergo forced decay, producing more delayed neutrons. It is importantto note that the precursor is not being transported in the simulation and only affectsthe simulation through the delayed neutrons that emits. To begin a transient Monte Carlo, an initial transient source distribution will beconstructed using the converged source distribution from an eigenvalue calculation,3
Precursor N e u t r o n p o p u l a t i o n F o r c e d p r e c u r s o r d e c a y F o r c e d p r e c u r s o r d e c a y F o r c e d p r e c u r s o r d e c a y F o r c e d p r e c u r s o r d e c a y F o r c e d p r e c u r s o r d e c a y P r e c u r s o r c r e a t i o n Figure 3.2: Schematic representation of the forced decay scheme for the precursors,where the precursor is forced to decay at the beginning of each time bin, so thereare scored in every time interval.with k eff ∼
1. To assess the convergence of the source distribution in the criticalitycalculation, OpenMC code has the capability of define a suitable spatial mesh andmonitor the Shannon entropy. There are two methods to create an initial particlesource. The first method is to transform the converged neutron source into a mixsource comprised of neutrons and precursors. The second method consists of sam-pling the initial neutrons and precursors using appropriate tallies after the eigenvaluecalculation. For the first stage in the development of this work, when a 1-group andmonoenergetic system was studied, the first method is an acceptable choice. As itwas shown in Section 2.1.2 and according to Eq. (2.4), the precursor concentrationfor one group at stationary state given by ( ∂C/∂t = 0) is C ( r ) = β i λ i ν Σ f ψ ( r ) . (3.14)To determine the fraction of neutrons in position r the following relation is useful n ( r ) n ( r ) + C ( r ) = v ψ ( r ) v ψ ( r ) + β i λ i ν Σ f ψ ( r ) = 11 + β i λ i νv Σ f . (3.15)This relation is valid only for constant neutron energy. For the mono-energeticsystem studied in this work and using the parameters shown in Table C.1, it is4obtained that for every neutron there are about 10 precursors and that the fractionof prompt neutrons in steady state is 0 .
08 %.For the second method [3] the energy dependent initial source is sampledfrom an eigenvalue calculation. For the number of neutrons, the estimator used was n ( r , E ) = (cid:90) π ψ ( r , Ω , E ) v ( E ) d Ω , (3.16)while for the precursors the estimator utilized was C i, ( r ) = (cid:90) π (cid:90) ∞ β i ( r , E ) ν ( r , E )Σ f ( r , E ) λ i ψ ( r , Ω , E ) dEd Ω , (3.17)where ψ is the flux sampled by an already existing flux tally in OpenMC.It is important to mention that the probability distribution for a precursorcreated in a fission event (shown in Sec. 3.3.3) is different than the one for a precursorcreated from the steady state distribution. This is because the precursors haveundergone a portion of its decay before t = 0. The different precursors with differentdecay constants result in a steady state group distribution given by P i = λ b λ i β i β , for precursor group λ b λ i I i , for individual precursor , (3.18)where λ b is the inversely weighted decay constant defined as λ b = β (cid:80) i β i λ i , for precursor group1 (cid:80) i I i λ i , for individual precursor , (3.19)This difference in the probability distributions is implemented in the code accordingto the time of creation of the precursor. When using the “forced decay” method the precursors always survive after they decayinto delayed neutrons. This means that the population of precursors is continuously5increasing, so population control for precursors must be implemented. The methodimplemented in the OpenMC code in this work is the
Combing method [38], whichwas originally developed for stationary Monte Carlo simulations. The idea of thismethod is to preserve the total statistical weight while mantaining a fixed number ofparticles. In the context of this work, keeping constant the number of particles servesfor two purposes: i) variance reduction and ii) reduced computing time by keepingthe population size approximately constant. If the system is super-critical, combingprevents the unlimited growth of the population, while if the system is sub-critical,keeps the simulation running by preventing the population from dying. If the systemis critical, combing prevents the divergence of the population due to fluctuations offission chains [39].If there are K particles at the end of a time interval and the objective is tocomb them to M particles. These K particles will be combed into M using a combwith M teeth. Figure 3.3 shows an example situation with K = 6 and M = 4. TheFigure 3.3: Diagram of the application of the combing method for 4 particles oftotal weight W combed into M = 3. The particles kept by the comb are particle 1,particle 3 and particle 4, each with weight W/ W = K (cid:88) i =1 = w i . (3.20)The comb teeth are equally spaced with the position of the teeth randomly selectedas t m = ξ WM + ( m − WM . (3.21)Each time a tooth hits interval i , the i -th article is duplicated ad assigned a weight w (cid:48) i = WM , (3.22)6where w (cid:48) i is the weight after combing. Defining the integer j by j < w i W/M ≤ j + 1 , (3.23)it can be seen that either j or ( j + 1) teeth of a comb with a pitch of W/M will hitan interval of length w i . In particular, the probability of j teeth fall in an interval i is p i,j = j + 1 − w i MW , (3.24)while the probability that j + 1 teeth fall in interval i is p i,j +1 = w i MW − j. (3.25)The expected weight for a single particle after combing is w (cid:48) i = p i,j j WM + p i,j +1 ( j + 1) WM = w i W/M WM = w i , (3.26)this implies that the combing preserves the total weight because after combing eachparticle is asigned a weight w (cid:48) i = W/M and since there are M particles, the totalweight is preserved. In this work both the neutron and precursor populations arecombed separately and for the monitoring of the precursor population the timedprecursor weight (Eq. (3.8)) or the expected delayed neutron weight (Eq. (3.10)) canbe used. hapter 4Results and discussion The time dependence in neutron transportation, including β -delayed neutron emis-sion from fission products, added in this work to the original OpenMC code weretested and the results are discussed in this chapter. This modified version of thementioned code will be denoted as Time-Dependent OpenMC or OpenMC(TD).In Section 4.1, the inclusion of time dependence and individual precursorsin OpenMC(TD) was evaluated. Related to time dependence, the tests made were:i) time tally (see Sec. 3.2.3), by scoring time dependent quantities in a fixed sourcecalculation in a subcritical configuration for the RA-6 reactor (see Sec. 4.1.1), ii)time boundary Monte Carlo simulation (see Sec. 3.2.2 and Sec. 3.2.1), by transportingneutrons in a fixed source calculation and in a Monte Carlo simulation divided in timeintervals (see Sec. 4.1.2), and iii) scoring of time dependent quantities in a simulationdivided in time intervals (see Sec. 4.1.3). The inclusion of individual precursors leadto a discussion about the β -delayed neutron activity comparing the standard 6-groupprecursor structure and the 50 individual precursor structure studied in this work (seeSec. 4.1.4). Likewise, the β -delayed average neutron energy for the 8-group precursorstructure in OpenMC(TD) was compared with the neutron spectrum energy for theJEFF-3.1.1 8-group precursor structure (see Sec. 4.1.5).In Section 4.2 a monoenergetic fissile system was simulated considering1-group precursor structure. Three configurations were studied and discussed: sub-critical (See Sec. 4.2.1), critical (Sec. 4.2.2) and reactivity insertion (See Sec. 4.2.3).378In Section 4.3 an energy dependent system using U was simulated consid-ering different precursor structures. Two configurations were studied and discussed:subcritical (see Sec. 4.3.1) and supercritical (see Sec. 4.3.2).In Section 4.4 an energy dependent and light-water moderated system using
U was simulated using different precursor structures and criticality configurations.Afterwards, the 6-group precursor structure effective multiplication factor was com-pared to the 50 individual precursor structure effective multiplication factor (seeSec. 4.4.1). Finally, the following cases were studied and discussed: i) comparisonbetween 6-group and 50 individual precursor structure in a critical configuration(see Sec. 4.4.2), and ii) comparison between 6-group, 50 individual and 40 individualprecursor structure in a critical configuration (see Sec. 4.4.3).Simulations were run at CSICCIAN (spanish acronym for Simulation andCalculation Center in Nuclear Sciences and Applications) clusters from the ChileanNuclear Energy Commission, its specifications are shown in Appendix D, along witha summary of the simulations presented in this work.
As explained in Sec. 3.2, there were some starting points that needed to be addressedin order to include time dependency in a Monte Carlo transport code. In short: i)time is explicitly added by means of time label to the particles, ii) the total simulationtime is divided in discrete time intervals and, iii) a new filter is added, so the codehas the capability to score time-dependent quantities. In order to check the correctimplementation of these characteristics into the code, three tests were conductedprior to the inclusion of the precursors and delayed neutrons.At the time of the writing of this thesis, measured β -delayed neutron energyspectra in databases [1, 10] were available only for 34 precursors [14]. In this workthe average energy of the β -delayed neutron was used for each individual precursor(see Sec. 4.1.5).9 Time [ms] N eu t r on f l u x [ a r b . un .] OpenMCMCNP6EXP
Figure 4.1: Time evolution of the neutron flux in a subcritical configuration for theRA-6 reactor, obtained from running a fixed source simulation in both OpenMC(TD)and MCNP. Results were scored for 10 ms. Both codes are in good agreement withthe experimental benchmark result.The tallying capabilities of OpenMC were expanded and a time filter wasadded to monitor the time evolution of any of the tallies already present in the code.In order to examine the proper functioning of this filter, MCNP and OpenMC(TD)were used to estimate the time evolution of the neutron flux in the RA-6 reactorusing the pulsed method [40]. In this method, a burst of neutrons is injected into asubcritical system and then the decay of the prompt neutron flux as a function oftime is observed. Since the phenomena being studied is the prompt neutron decay thecontribution from delayed neutrons can be neglected from point kinetics Eq. (2.17),which in that case reads, ∂∂t s ( t ) = ρ − β eff Λ eff s ( t ) . (4.1)The solution to Eq. (4.1) is given by s ( t ) = s e α t , with α = ρ − β eff Λ eff (4.2)0where s is the initial flux density and the decay constant α is the decay constant ofthe neutron population.MCNP and OpenMC(TD) were used to simulate the neutron source andthen the prompt neutron decay was scored during 10 ms. The flux as a func-tion of time obtained with both codes was then compared with the experimentalresults for the decay constant from the graphite reflected RA-6 benchmark from[33, 34]. Fig. 4.1 shows results obtained, where blue (red) crosses (x marks) denoteOpenMC(TD) (MCNP) results. Dashed curve denote the benchmark value. It canbe observed good agreement between the decay constants from fit parameters andthe result from the benchmark. OpenMC MCNP Benchmark α [s − ] − − − time filter imple-mented works as expected and OpenMC(TD) can score time dependent quantitiesin fixed source calculations. Another modification needed to be the implemented in the code is the division of thetotal simulation time in discrete time intervals. It is important to check that thereare no errors in the crossing of time intervals. To do this, neutron transport in themonoenergetic fissile system described in Appendix C was studied. Since the purposeof this test was only to check for errors in the particle transport when dividing thesimulation in time intervals, fission reactions were not considered. Results obtainedare shown in Fig. 4.2, where it can be seen that the neutron flux obtained whenthe simulation is divided in discrete time intervals is the same when a regular fixed1 - · Time [s] - F l u x [ a r b i t r a r y un i t s ] Time intervalsContinuous
Figure 4.2: Neutron flux as a function of time in a simple transport problem. In red,the neutron flux obtained for a non-transient fixed source simulation is shown. Inblue the neutron flux obtained from a transient simulation divided in time intervalsis shown. Both results are equivalent.source calculation is performed, thus, the transport logic is correct and OpenMC(TD)correctly transports neutrons across time intervals.
Test 3 was a combination of tests 1 and 2, i.e., flux scoring as a function of timein a subcritical configuration when the simulation was divided in time intervals.The transport problem studied was the monoenergetic fissile system detailed in Ap-pendix C. The advantage of studying a system like this one is that the time evolutionof the neutron flux can be described by an analytical expression, making direct itsvalidation.First, the configuration was made subcritical by increasing the absorptioncross section to Σ a = 0 . − , while mantaining the total cross section constant.A criticality calculation with 10 neutrons, 5000 batches and 300 skipped cycles for2 - - F l u x [ a r b i t r a r y un i t s ] MCAnalytical
Figure 4.3: Time evolution of the neutron flux for the monoenergetic system studied.Results obtained using OpenMC(TD) are shown in blue, while the fit to the pointkinetics solution given by Eq. (4.2) is shown in red.this configuration gives k eff = 0 . ± . neutrons for 300 ms using a time interval of 1 ms and 10 batches. Fig. 4.3shows a comparison between prompt neutron flux obtained from the Monte Carlosimulation and the analytical solution obtained using Eq. (4.2) and the parametersfor this system given in Table C.1. The fitted time constant parameter obtained forthe decay of the prompt neutron flux is α fitted = − . − , while the calculatedvalue is given by α teo = − . − . A summary of the obtained results is shown inTable 4.2. Both values are in excellent agreement with each other (1 .
2% difference).Therefore, the scoring of time dependent quantities when the simulation is dividedin time intervals works correctly.
Calculated Fitted ∆ Decay constant Decay constant α [s − ] − . − . − . The purpose behind the activity calculation for individual precursors and its compar-ison to the 6-group activity was to verify the suitability of the 50 individual precur-sors chosen for the emission of the β -delayed neutrons as part of the new capabilitiesOpenMC(TD) code. This test was necessary because if there were differences inresults obtained for the time evolution of the neutron flux using the N -group struc-ture or in the individual precursors, it was relevant to know if the activity of the β -delayed neutron emission was the cause of these eventual discrepancies.The calculated activity for the 6-precursor groups, denoted by A ( t ), isgiven by A ( t ) = (cid:88) i =1 a i exp( − λ i t ) , (4.3)where a i = β i /β is the i -th group relative abundance and λ i is the i -th group decayconstant (see Table E.1). Conversely, the calculated activity for the 50 individualprecursors, denoted by A ( t ) reads A ( t ) = (cid:88) i =1 I i exp( − λ i t ) , with I i = CY i P n,i ν d , (4.4)where I i is the i -th precursor importance as defined in Eq. (3.3) (see Table E.2) and λ i is the i -th precursor decay constant (See Sec. 2.3.1). The calculated activity for 6precursor groups and 50 individual precursors is shown in Fig. 4.4. In blue, A ( t ) isshown, while A ( t ) is shown in red. As it can be seen, both activities are equivalent.Quantitatively, comparing β -delayed neutron emission for A ( t ) and A ( t ) up to100 s, it is obtained that (cid:82) A ( t ) dt (cid:82) A ( t ) dt = 0 . . (4.5)This indicates that adding the remaining 219 precursors only contributes to 0 . - -
10 1 - de l a y ed neu t r on a c t i v i t y b Figure 4.4: β -delayed neutron activity for 6 precursor groups and 50 individualprecursors is shown. In blue, A ( t ) is shown, while A ( t ) is shown in red. As it canbe seen, both activities are equivalent. At the time of the writing of this work there exists experimental measurements foronly 34 β -delayed neutron energy spectra. This data was compiled and completedby Brady in 1989 [14]. The remaining β -delayed neutron energy spectra presentin ENDF/B-VIII.0 comes from QRPA calculations [1]. Given that the capabilitiesadded to the OpenMC code allows to run simulations using up to 269 individualprecursors and with the intention of having these precursors on the same footingregarding the β -delayed neutron emission energies, it was decided that the averageenergy for the delayed neutron emission would be used. Nevertheless, if the β -delayedneutron energy spectra databases were updated in the future, its inclusion could beeasily implemented in the code.Since the average energy for the β -delayed neutron emission was used, it wasimportant to verify that the results obtained for the time evolution of the neutron5flux when using the delayed neutron average energies were equivalent to samplingthe delayed neutron energy from the corresponding spectra. To this end, a transientsimulation using OpenMC(TD) in a subcritical configuration was run. Subcriticalitywas achieved by decreasing the U density from δ U = 3 . × − (atoms/bcm) to δ U = 3 . × − (atoms/b cm), while mantaining the dimensions of thebox constant, obtaining an effective multiplication system of k eff = 0 . ± . µ s each. Population control wasapplied at the end of each interval.Results obtained from transient Monte Carlo simulation using OpenMC(TD)for the time evolution of the neutron flux when the delayed neutron energy was sam-pled from spectra are shown in red in Fig. 4.5, while results obtained when using theaverage energy for the delayed neutron emission are shown in blue. From Fig. 4.5it can be seen that the time evolution of the neutron flux obtained with transientMonte Carlo code OpenMC(TD) using the average delayed neutron energy and de-layed energy sampled from spectra are equivalent. Once the preliminary work described in Section 4.1 was completed, the new capa-bilities added to the OpenMC code, namely division of the simulation in discretetime intervals, scoring of time dependent quantities, forced decay of precursors andpopulation control, were tested in the monoenergetic fissile system described in Ap-pendix C. The objective of this section is to lay the groundwork for the study of thedelayed neutron emission from individual precursors for transient calculations usingOpenMC.Prior to transient simulations with Monte Carlo code OpenMC(TD) pre-sented in this section, a non-transient standard steady state criticality calculationwas done with 10 neutrons, 3000 batches and 200 skipped cycles, using OpenMC.The effective multiplication factor obtained was k eff = 1 . ± . - -
10 1 N o r m a li z ed neu t r on f l u x EnergyMean energy
Figure 4.5: Time evolution of neutron flux for a subcritical configuration with k eff =0 . ± . µ s each. Results of neutron flux whendelayed neutron energy was sampled from spectra are shown in red (behind the bluecurve), while results obtained when using the average energy for delayed neutronemission are shown in blue. It can be seen that both results are equivalent.wards, the initial transient source was created as described in Sec. 3.4, with 10 neutrons and 9 × precursors. This initial transient source was used in subcriti-cal (See Sec. 4.2.1), critical (See Sec. 4.2.2) and reactivity insertion (See Sec. 4.2.3)configurations, presented in the following subsections, in order to start the transientsimulation.In this section, the code input was the macroscopic absorption cross section,Σ a , which was suitably modified in order to produce reactivity changes in the mo-noenergetic fissile system (critical, subcritical or supercritical configurations). Thus,the output values (observables) were the effective multiplication factor k eff and thetime evolution of the neutron flux φ ( t ), which were compared with point kineticscalculations.7 The code was firstly tested in a subcritical configuration. Subcriticality was achievedby increasing the absorption cross section from Σ a = 0 . − to Σ a = 0 . − .Total cross section Σ t was kept constant, then the effective multiplication of thesystem k eff = 0 . ± . a is equivalent to decrease the density of the fissile material δ f of the system.The simulation was run using 60 batches and the total simulation time was50 s, divided in 500 time intervals of 100 ms each one. At the end of each time intervalpopulation control was applied, using the technique explained in Sec. 3.5. Resultsobtained from transient Monte Carlo simulation using OpenMC(TD) are shown inblue in Fig. 4.6, meanwhile the point kinetic solution of the neutron population as afunction of time is shown in red. From Fig. 4.6 it can be seen that the time evolutionof the neutron population calculated using transient Monte Carlo code OpenMC(TD)and point kinetics solution using Eq. (B.4) are equivalent.Quantitatively, from Fig. 4.6 the reactivity value ρ can be obtained asa fitted parameter of Eq. (B.2). This ρ fit = − . ρ = ( k eff − /k eff = − . . This new time dependent capability addedto OpenMC allows the observation of the slow decay of the neutron population dueto the β -delayed neutron emission.The reactivity value is usually obtained by running a criticality Monte Carlocalculation. In this work, using OpenMC(TD), this value can be obtained by fittingEq. (B.4) to the time evolution of the neutron population. A summary of the resultsobtained is shown in Table 4.3. See for instance Fig. 4.3 where in a non-transient standard Monte Carlo fixed source calculation,using both MCNP or OpenMC, the neutron population extinguishes in ∼
50 ms. Calculated Fittedreactivity reactivity ρ [pcm] − − - -
10 1 N o r m a li z ed neu t r on popu l a t i on OpenMC(TD)Point kinetics
Figure 4.6: Time evolution of the neutron population for a monoenergetic system ina subcritical configuration ( k eff = 0 . ± . In the second test, transient analysis was done in a critical configuration, with k eff =1 . ± . N o r m a li z ed neu t r on popu l a t i on OpenMC(TD)Point kinetics
Figure 4.7: Time evolution of the neutron population for a monoenergetic system in acritical configuration ( k eff =1 . ± . ρ can be obtained as afitted parameter of Eq. (B.2). This ρ fit = 0 . ρ = ( k eff − //k eff = 0 . Calculated Fittedreactivity reactivity ρ [pcm] 10(3) 13(70)Table 4.4: Results obtained for the reactivity of the monoenergetic simulated systemin a critical configuration using 1-group precursor structure. · N eu t r on popu l a t i on [ a r b . un i t s ] OpenMC(TD)Point kinetics
Figure 4.8: Time evolution of the neutron population for a monoenergetic systemobtained using OpenMC(TD). The system is initially in a critical configuration, then,at t = 10 s a reactivity of 211 pcm is inserted. After 30 s the system is brought backto critical configuration.The last case studied was a mixture of the two cases previously presented.First, the configuration was critical, then reactivity is inserted and afterwards, theconfiguration is brought back to critical by a negative reactivity insertion. Thissystem is a first approximation to simulate the operation of a nuclear reactor. Con-cretely, these reactivity insertions were simulated by changing the absorption crosssection while keeping the total cross section constant.1In Fig. 4.8 the reactivity insertion case is shown. For the first 10 s theconfiguration was kept critical with Σ a = 0 . − . Then, for 10 s < t <
40 s theabsorption cross section was reduced from Σ a = 0 . − to Σ a = 0 . − ,inserting a positive reactivity of 211 pcm, thus making the configuration supercritical.This fast increase in neutron population is known as prompt jump . In t = 40 s thesystem is brought back to Σ a = 0 . − . The neutron population stops growingand decreases rapidly. This fast change in neutron population is known as promptdrop . In t = 40 s the neutron population is almost three times the initial neutronpopulation. The final state of the configuration is slightly supercritical.The simulation was run for 25 batches and the simulation time was dividedin 5000 time intervals of 10 ms each and population control was applied at the endof every interval. Results obtained are shown in Fig. 4.8, where the time evolutionof the neutron population calculated from point kinetic equations is also shown.From Fig. 4.8 it can be seen that the time evolution of the neutron populationcalculated using transient Monte Carlo simulation using OpenMC(TD) and pointkinetics solution using Eq. (B.4) are equivalent.It is important to notice that in the reactivity insertion case, both promptjump and prompt drop can be studied in detail given that short time intervals of10 ms were used in the simulation. This new Monte Carlo capabilitity, implementedin this work, allows to reduce time windows as much as desired, so parameters asthe Rossi- α [41] can be calculated. U system
After the new capabilities added to the code were successfully tested for the monoen-ergetic system described in the previous Section 4.2, the following study involvedtesting the code in a system with continuous, energy-dependent cross sections (i.e.not monoenergetic). The objective of this section is to simulate a more realisticsystem, but at the same time keeping it simple enough to compare to the pointkinetics model, whenever is possible. In order to do this, the material of the boxfrom the preceding section was made of pure
U, using the continuous energy crosssections from JEFF-3.1.1 [10] nuclear data library and the geometry was surrounded2by vacuum.Prior to the transient simulations with Monte Carlo code OpenMC(TD)presented in this section, a non-transient standard steady state criticality calculationwas done with 10 neutrons, 5000 batches and 300 skipped cycles, using OpenMC.The effective multiplication factor obtained was k eff = 1 . ± . neutrons and 9 × precursors. This initial transient source was used in subcriticaland supercritical tests, presented in the following subsections, in order to start thetransient simulation. Intentionally, the critical configuration was not considered inthis set of tests because the main objective of this part of the work was to examinewhether the code had the capability to resolve fast changes in the neutron flux.In this section, the code input was the density of the fissile material, whichwill be denoted as δ U . Since the box made of pure U, the only two ways to insertreactivity to the system are: i) by changing the box dimensions or, ii) by changingthe density of the fissile material. The latter method was chosen and the dimensionsof the box were kept constant throughout the different cases. Thus, the output value(observable) was the effective multiplication factor k eff and the time evolution of theneutron flux φ ( t ), like in the previous section. Since this is not a monoenergeticsystem, Eq. (C.1) for the effective generation time no longer holds. In consequence,for the calculation of Λ and β eff , a simulation of the system in MCNP was made,given that this code can estimate these parameters using the weighted adjoint flux.These two quantities were then compared with the fitted parameters from Eq (B.4),which is the solution to the point kinetics equations.Different group structures were simulated in this section. When it waspossible, the energy of the β -delayed neutrons was taken from a distribution (JEFF-3.1.1). Otherwise, the average energy was used for each precursor or group (ENDF-B/VIII.0). For comparison purposes, a simulation using the energy distribution andthe average energy from the first group were also studied.3 The first case studied was a subcritical configuration. The system was made sub-critical by decreasing the
U density from δ U = 4 . × − (atoms/b cm) [3] to δ U = 4 . × − (atoms/b cm), while mantaining the dimensions of the box con-stant, making the effective multiplication factor of the system k eff = 0 . ± . i) First group with energy distribution from JEFF-3.1.1 The first precursor group, characterized by a half-life T / = 55 . β -delayed neutron energy spectrum from JEFF-3.1.1. is shown in Fig. 4.9. In Appendix A, the β -delayed neutron energy spectra forall 8 groups can be found. × Energy (eV) − × d N / d E ( / e V ) xx Mean = 212300.9 eV
Figure 4.9: Group 1 β -delayed neutron energy spectrum from JEFF-3.1.1.This simulation was run using 22 batches and the total simulation time was0 . µ s, and then for t > µ s the decay of the neutron flux stabilizes. Fit parameters β ± Λ − ± − − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x x x Figure 4.10: Study i). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, and then the decay of the neutron population slows.Inset figure shows the prompt drop zoomed for the first 5 µ s.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 5 . ∼ .
1% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 5 . . . β eff [pcm] 644(6) 648(38) 1%Table 4.5: Study i). Values of the parameters obtained from running anOpenMC(TD) transient simulation in a subcritical configuration. The calculatedvalues for the parameters were calculated using MCNP, while the OpenMC(TD)parameters were obtained by fitting Eq. (B.4) to the time evolution of the neutronflux. ii) First group with average energy from JEFF-3.1.1 The first precursor group was simulated, but in this case the delayed neutron wasemitted with the average energy of the first group energy distribution reported fromJEFF-3.1.1. This energy is ¯ E g = 212 .
31 keV and it was calculated as the weightedaverage per eV from the distribution shown in Fig. 4.9.This simulation was run using 3 batches and the total simulation time was0 .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.11, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.11 the prompt drop can be seenfor the first 5 µ s, and then for t > µ s the decay of the neutron flux stabilizes.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 5 . ∼ .
1% difference between both quantities. The fitted effective delayedneutron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Fit parameters β ± Λ − ± − − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x x Figure 4.11: Study ii). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, and then the decay of the neutron population slows.Inset figure shows the prompt drop zoomed for the first 5 µ s. Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 5 . . . β eff [pcm] 644(6) 666(34) 3 . iii) -group with average energy from ENDF-B/VIII.0 A 1-group precursor structure was simulated. Delayed neutrons were emitted withthe average energy of the 6-group precursor structure from ENDF-B/VIII.0. This7energy is ¯ E g = 501 .
31 keV and was calculated as the weighted average of thereported average energies per group, according to¯ E g = (cid:88) i =1 β i β ¯ E i , (4.6)where E i is the average energy for i -th group, see Table 2.1.This simulation was run using 3 batches and the total simulation time was0 .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.12, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.12 the prompt drop can be seenfor the first 5 µ s, and then for t > µ s the decay of the neutron flux stabilizes.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 5 . ∼ .
7% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 5 . . . β eff [pcm] 644(6) 602(36) 6 . Fit parameters β ± Λ − ± − − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x x Figure 4.12: Study iii). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, and then the decay of the neutron population slows.Inset figure shows the prompt drop zoomed for the first 5 µ s. iv) -group with energy distribution from JEFF-3.1.1 An 8-group precursor structure was simulated. Delayed neutrons energies were ran-domly sampled from one of the energy distributions from JEFF-3.1.1, shown inAppendix A.This simulation was run using 3 batches and the total simulation time was0 .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.13, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.13 the prompt drop can be seen9for the first 5 µ s, and then for t > µ s the decay of the neutron flux stabilizes. Fit parameters β ± Λ − ± − − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x x Figure 4.13: Study iv). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, and then the decay of the neutron population slows.Inset figure shows the prompt drop zoomed for the first 5 µ s.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 5 . ∼ .
1% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . v) -group with average energy from ENDF-B/VIII.0 A 6-group structure was simulated. Delayed neutrons energies were randomly sam-pled according to β i /β , from the listed average energies of the six precursor groups0 Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 5 . . . β eff [pcm] 644(6) 660(60) 2 . .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.14, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.14 the prompt drop can be seenfor the first 5 µ s, and then for t > µ s the decay of the neutron flux stabilizes.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 5 . ∼
1% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . vi) individual precursors with average energies from ENDF-B/VIII.0 A 50 individual precursor structure was simulated. Delayed neutrons were randomlysampled according to its importances I i =( CY i P n,i ) /ν d , from the calculated averageenergies of the 50 individual precursors used in this work (see Table E.2).This simulation was run using 3 batches and the total simulation time was1 Fit parameters β ± Λ − ± − − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x x Figure 4.14: Study v). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, and then the decay of the neutron population slows.Inset figure shows the prompt drop zoomed for the first 5 µ s. Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 5 . . β eff [pcm] 644(6) 602(57) 6 . .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation using2OpenMC(TD) are shown in blue in Fig. 4.15, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.15 the prompt drop can be seenfor the first 5 µ s, and then for t > µ s the decay of the neutron flux stabilizes. Fit parameters β ± Λ − ± − − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]0.20.40.60.811.21.4 N o r m a li z ed neu t r on f l u x x Figure 4.15: Study vi). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, and then the decay of the neutron population slows.Inset figure shows the prompt drop zoomed for the first 5 µ s.3Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 5 . ∼ .
3% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 5 . . . β eff [pcm] 644(6) 602(57) 6 . Now a supercritical configuration was studied. The system was made supercriticalby increasing the
U density from δ U = 4 . × − (atoms/b cm) to δ U =4 . × − (atoms/b cm), while mantaining the dimensions of the box constant,making the effective multiplication factor of the system k eff = 1 . ± . i) First group with energy distribution from JEFF-3.1.1 The first precursor group, characterized by a half-life T / = 55 . β -delayed neutron energy spectrum from JEFF-3.1.1. is shown in Fig. 4.9, and the remaining delayed neutron group spectra can befound in Appendix A.The simulation was run using 10 batches and the total simulation time was0 . µ s, and then for t > µ s the growth of the neutron flux stabilizes.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 6 . ∼ .
2% difference between both quantities. The fitted effective delayed neutronfraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff =0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 6 . . . β eff [pcm] 651(6) 666(56) < Fit parameters β ± Λ − ± − − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x x Figure 4.16: Study i). Time evolution of the neutron flux in the studied supercriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD) the configu-ration is made supercritical by increasing δ U . The time evolution of the neutronflux is shown in blue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutron populationslows. Inset figure shows the prompt jump zoomed for the first 10 µ s. ii) First group with average energy from JEFF-3.1.1 The first precursor group was simulated, but in this case the delayed neutron wasemitted with the average energy of the first group energy distribution reported fromJEFF-3.1.1. This energy is ¯ E g = 212 .
31 keV and it was calculated as the weightedaverage per eV from the distribution shown in Fig. 4.11.This simulation was run using 3 batches and the total simulation time was0 .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.17, while the fit obtained by adjusting6the results to Eq (B.4) are shown in red. In Fig. 4.17 the prompt jump can be seenfor the first 1 µ s, and then for t > µ s the growth of the neutron flux stabilizes. Fit parameters β ± Λ − ± − − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x x Figure 4.17: Study ii). Time evolution of the neutron flux in the studied supercriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD) the configu-ration is made supercritical by increasing δ U . The time evolution of the neutronflux is shown in blue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutron populationslows. Inset figure shows the prompt jump zoomed for the first 10 µ s.7Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 6 . ∼ .
2% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 6 . . . β eff [pcm] 651(6) 666(63) 2 . iii) -group with average energy from ENDF-B/VIII.0 A 1-group precursor structure was simulated. Delayed neutrons were emitted withthe average energy of the 6-group precursor structure from ENDF-B/VIII.0. Thisenergy is ¯ E g = 501 .
31 keV and it was calculated as the weighted average of thereported average energies per group, according to¯ E g = (cid:88) i =1 β i β ¯ E i , (4.7)where E i is the average energy for i -th group, see Table 2.1.This simulation was run using 3 batches and the total simulation time was0 .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.18, while the fit obtained by adjusting8
Fit parameters β ± Λ − ± − − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x x Figure 4.18: Study iii). Time evolution of the neutron flux in the studied supercriti-cal configuration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD) the configu-ration is made supercritical by increasing δ U . The time evolution of the neutronflux is shown in blue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutron populationslows. Inset figure shows the prompt jump zoomed for the first 10 µ s.the results to Eq (B.4) are shown in red. In Fig. 4.18 the prompt jump can be seenfor the first 1 µ s, and then for t > µ s the growth of the neutron flux stabilizes.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 6 . ∼ .
2% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 6 . . . β eff [pcm] 651(6) 637(35) 2 . iv) -group with energy distribution from JEFF-3.1.1 An 8-group precursor structure was simulated. Delayed neutrons energies were ran-domly sampled from one of the energy distributions from JEFF-3.1.1, shown inAppendix A.This simulation was run using 3 batches and the total simulation time was0 .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.19, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.19 the prompt jump can be seenfor the first 1 µ s, and then for t > µ s the growth of the neutron flux stabilizes.0 Fit parameters β ± Λ − ± − − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x x Figure 4.19: Study iv). Time evolution of the neutron flux in the studied supercriti-cal configuration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD) the configu-ration is made supercritical by increasing δ U . The time evolution of the neutronflux is shown in blue, while the fit obtained is shown in red. Between 0 < t < µ sthe prompt jump can be observed, an then the growth of the neutron populationslows. Inset figure shows the prompt jump zoomed for the first 10 µ s.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 6 . MCNP = 6 . ∼ <
1% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 6 . . < β eff [pcm] 651(6) 665(35) 2 . v) -group with average energy from ENDF-B/VIII.0 A 6-group structure was simulated. Delayed neutrons energies were randomly sam-pled according to β i /β , from the listed average energies of the six precursor groups(see Table 2.1). This simulation was run using 3 batches and the total simulationtime was 0 .
05 ms divided in 500 time intervals of 100 ns each. Population control(see Sec. 3.5) was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.20, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.20 the prompt jump can be seenfor the first 1 µ s, and then for t > µ s the growth of the neutron flux stabilizes.2 Fit parameters β ± Λ − ± − − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x x Figure 4.20: Study v). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, and then the decay of the neutron population slows.Inset figure shows the prompt drop zoomed for the first 10 µ s.Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 6 . ∼ .
2% difference between both quantities. The fitted effective delayed neu-tron fraction obtained was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 6 . . . β eff [pcm] 651(6) 635(38) 2 . vi) individual precursors with average energies from ENDF-B/VIII.0 A 50 individual precursor structure was simulated. Delayed neutrons were randomlysampled according to its importances I i = ( CY i P n,i ) /ν d , from the calculated averageenergies of the 50 individual precursors used in this work (see Table E.2).This simulation was run using 3 batches and the total simulation time was0 .
05 ms divided in 500 time intervals of 100 ns each. Population control (see Sec. 3.5)was applied at the end of each interval.Results obtained from the transient Monte Carlo simulation usingOpenMC(TD) are shown in blue in Fig. 4.20, while the fit obtained by adjustingthe results to Eq (B.4) are shown in red. In Fig. 4.21 the prompt jump can be seenfor the first 1 µ s, and then for t > µ s the growth of the neutron flux stabilizes. Fit parameters β ± Λ − ± − − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x OpenMC(TD)Fit x − × Time [s]00.20.40.60.811.21.41.61.82 N o r m a li z ed neu t r on f l u x x Figure 4.21: Study vi). Time evolution of the neutron flux in the studied subcriticalconfiguration. The initial transient source is prepared in a critical state and at thebeginning of the transient Monte Carlo simulation using OpenMC(TD), the systemis made subcritical by decreasing δ U . The time evolution of the neutron flux isshown in blue, while the fit obtained is shown in red. Between 0 < t < µ s theprompt drop can be observed, after which the decay of the neutron flux slows. Insetfigure shows the prompt drop zoomed for the first 10 µ s.5Quantitatively, the fitted effective neutron generation time obtained wasΛ fitted = 5 . MCNP = 6 . ∼ .
2% difference between both quantities. The fitted effective delayed neu-tron fraction was β ( fitted ) eff = 0 . β ( MCNP ) eff = 0 . Parameter Calculated Fitted ∆ Unit MCNP OpenMC
Λ [ns] 6 . . . β eff [pcm] 651(6) 621(36) 4 . . µ s are used in the simulation. This new Monte Carlo capability was shownto work in a monoenergetic system and here it is shown that also works in a varyingenergy system. Lastly, the OpenMC(TD) code also correctly describes the behaviourof the system after the prompt jump, when the neutron flux changes slowly.In this section, values for: i) nuclear cross sections, ii) mean energies, iii)cumulative yields, iv) probability of delayed neutron emission and, v) decay constantswere taken from nuclear databases JEFF-3.1.1 and ENDF-B/VIII.0. In this regard,JEFF-3.1.1 has reported the neutron energy spectra for each of the eight groups, butit does not have the delayed neutron energy spectra for the 269 individual precursors.As it was seen in Section 2.3.2, there are discrepancies between both databases.In this work, it was necessary to use the cumulative yields from JEFF-3.1.1, andthe probability of delayed neutron emission and average delayed neutron energyfrom ENDF-B/VIII.0. These discrepancies are finally reflected in the value obtainedfor the effective delayed neutron fraction ( ¯ β eff = 658(26) pcm for JEFF-3.1.1 and¯ β eff = 602(29) pcm for ENDF-B/VIII.0, for the subcritical configuration; ¯ β eff =666(34) pcm for JEFF-3.1.1 and ¯ β eff = 631(21) pcm for ENDF-B/VIII.0, for thesupercritical configuration). This shows how important is that every database countswith good and better nuclear data for individual precursors, either average energiesor energy spectra. In this section the fast system studied in Section 4.3 was modified by including aneutron moderator surrounding the
U. The β -delayed neutron emission now wasproduced by individual precursors and results obtained were compared when emissionis from the 6-group precursor structure.Comparisons were made between simulations using the 6-group structureand 50 individual precursors, such as i) effective multiplication factor for a critical7system (see Sec. 4.4.1), and ii) time evolution of the neutron flux in a transientsimulation (see Sec. 4.4.2).As a final test the 10 most important precursors were removed from the50 individual precursors, in order to account for its effect on the time evolution ofthe neutron flux in comparison with the 6-group structure (see Sec. 4.4.3).The configuration was surrounded with a 4 .
29 cm-thickness water modera-tor and made critical by setting the
U density to δ U = 3 . × − (atoms/bcm). The continuous energy cross sections used were from JEFF-3.1.1 [10]. Thedimensions of the box remained the same. Prior to the transient simulations pre-sented in this section, a non-transient standard criticality calculation was run with10 neutrons, 5000 batches and 300 skipped cycles using OpenMC. The effectivemultiplication factor obtained was k eff = 1 . ± . δ U , the delayed neutron energy (sampled or averaged fromspectra) and the number of precursors used (50 or 40). Reactivity was inserted usingthe same method described in Sec. 4.3. The output values (observables) were theeffective multiplication factor k eff and the time evolution of the neutron flux φ ( t ), likein the previous section. Since this is no longer a 1-group precursor problem, there areno analytical solutions to the point kinetics equations. Nevertheless, resorting to thepoint kinetics approximation, a good estimation to the asymptotic decay constantfor the neutron flux [6] can be found using the equation α D = ¯ λρβ eff − ρ , (4.8)where α D is the asymptotic decay constant, ¯ λ is the average β -weighted decay con-stant , β eff is the effective delayed neutron fraction and ρ is the system reactivity.Regarding the choice of the effective delayed neutron fraction, the average delayedneutron yield obtained when using the data from JEFF-3.1.1 in Eq. (2.28) is ν d =1 . × − , while the value obtained when using ENDF/B-VIII.0 is ν d = 1 . × − .If the average neutron yield is taken to be ν = 2 . β = ν d /ν , ranges from β = 607 pcm to β = 780 pcm. In view of this, the Importance was defined in Eq. (3.3), see Sec. 3.3.2. The average β -weighted decay constant is given by ¯ λ = (cid:80) β i β λ i β eff = 700 pcm.The decay constant was ¯ λ = 0 . − . The reactivity of the system was obtainedas fitted parameter, and then compared to the reactivity obtained from the initialnon-transient criticality calculation ( ρ = ( k eff − /k eff ). As it was shown in Section 4.2 and Section 4.3, prior to every transient simulationwith Monte Carlo code OpenMC(TD), a non-transient, standard steady state criti-cality calculation with OpenMC must be done in order to create the initial transientsource, and assess the reactivity of the system. Since during the writing of thisthesis there are no codes able to perform a criticality calculation using individualprecursors as the source of β -delayed neutrons, in this work the capability to runcriticality calculations using individual precursors instead of the N -group structurewas also added to the OpenMC(TD) code. Criticality was achieved by mantaining U density at δ U = 3 . × − (atoms/b cm) obtaining k eff = 1 . k eff = 1 . precursors Difference k eff . . U cube when is ther-malized by surrounding it with a water moderator of a 4 .
29 cm thickness. Resultsfor 6-group structure and 50 individual precursors.The effective multiplication factor obtained for this system shows that thisconfiguration is slightly supercritical and both results are in good agreement witheach other, with a difference of 7 pcm among them.9
A transient simulation using OpenMC(TD) was ran for the previous system (seeSec. 4.4.1) in a critical configuration comparing the time evolution of the neutronflux obtained when 6-group and 50 individual precursor structures were used. Bothsimulations were run using 2 batches. Total simulation time was 4 s divided in400 time intervals of 10 ms each. Population control was applied at the end ofeach interval. The wall-clock time for the 6-group precursor simulation was about260 .
05 h, while for the 50 individual precursor simulation was 410 .
76 h.Results obtained from transient Monte Carlo simulation using OpenMC(TD)for the time evolution of the neutron flux, for the 6-group structure, are shown inblue in Fig. 4.22, while results obtained for the 50 individual precursor structure areshown in red. From Fig. 4.22 it can be seen qualitatively that both results show aslightly supercritical system, where the neutron flux increases slowly in time, whichis consistent with the effective multiplication factor of a near critical configuration.0 N o r m a li z ed neu t r on f l u x Figure 4.22: Time evolution of neutron flux in a water moderated box made of pure
U, simulated with OpenMC. Results for 6 groups are shown in blue, and results for50 individual precursors are shown in red. Both results show a slightly supercriticalsystem, where neutron flux increases slowly in time, consistent with the k eff of a nearcritical configuration.Now, analyzing Fig. 4.22, from a quantitative point of view, the reactivityvalue ρ can be obtained as fitted parameter of φ ( t ) ∼ e − α D t , for both 6-group and 50individual precursor structure . These ρ (6) fit = 0 . ρ (50) fit = 0 . ρ (6) = ( k (6) eff − /k (6) eff = 0 . ρ (50) = ( k (50) eff − /k (50) eff =0 . φ ( t ) ∼ e − α D t to the time evolution of the neutron population. A summary ofthe results obtained can be seen in Table 4.18.From examining the results obtained for the fitted parameters, it can benoticed that even when they are in good agreement with the calculated values, theypossess quite large uncertainties. This is due to the fact that neutron population for k eff > ρ >
0) shows an exponential growing behaviour. Simulations ran for 4 swhich was not enough time to reveal the exponential growing. To reduce the fitted The asymptotic decay constant α D was defined in Eq. (4.8). -group individualstructure structure ρ [pcm] 25(3) 32(3) ρ fit [pcm] 17(368) 35(347)Table 4.18: Results obtained for the reactivity of the water moderated energy de-pendent simulated system in a critical configuration using 6-group and 50 individualprecursor structure.uncertainties, the simulation time should increase to tens of seconds. This wouldincrease the wall-clock time of the simulation. For instance, a simulation time of50 s, would take 216 days (7 months and 6 days), beyond reach for the purposes ofthis thesis.Nevertheless, to further explore these uncertainties issues related to thesimulation time, the fitted reactivities from the subcritical (where the neutron pop-ulation was scored for 50 s, see Sec. 4.2.1) and critical (where neutron populationwas scored for 25 s, see Sec. 4.2.2) configurations of the monoenergetic fissile systemfrom Sec. 4.2, were obtained by taking into account the time evolution of neutronpopulation only for the first 4 s.For the monoenergetic system with a subcritical configuration (see Sec. 4.2.1),the previously fitted reactivity was ρ (50 s ) fit = − . ρ = − . ρ (4 s ) fit = − . ρ (50 s ) fit and ρ are in excellent agreement with each other, ρ (4 s ) fit shows a difference of2 .
85% with ρ . But the uncertainty of the calculated reactivity is 0 . ρ (4 s ) fit is almost 25 times its value.For the monoenergetic system with critical configuration (see Sec. 4.2.2),the previously fitted reactivity was ρ (25 s ) fit = 0 . ρ = 0 . ρ (4 s ) fit = 0 . ρ (4 s ) fit is almost 66 times its value.It is important to remark that OpenMC(TD) is stable for both configura-2tions and the rate of change for the neutron flux is consistent with the reactivitiesfrom the criticality calculations, considering that: i) continuous energy-dependentcross sections, ii) addition of a neutron moderator to the system and, iii) implemen-tation of individual precursors.In summary, the time evolution of neutron flux in a water moderated boxof U was obtained using 50 individual precursors and the results were consistentwith the initial calculated reactivities. The simulation was stable, and there was nodivergence of the neutron fission chains, which means that population control workedas expected.
The final study was a critical configuration, but in this case the 10 precursors withthe largest importances I i (see Eq. (3.3)) were removed from the previous individualprecursor structure. This means that a 40 individual precursor structure was usedfor this calculation. As in Sec. 4.4.2, the simulation was run using OpenMC(TD)for 2 batches and the total simulation time was 4 s, divided in 400 time intervalsof 10 ms each. Population control was applied at the end of each time step. Thewall-clock time for this simulation was about 319 .
65 h.Fig. 4.23 shows results obtained from transient Monte Carlo simulationusing OpenMC(TD) for the time evolution of the neutron flux. Results in blue arewhen 6 individual precursor structure was used, in red, when 50 individual precursorstructure was used, while in green when 40 individual precursor structure was used.In this case it can be seen that the time evolution of the neutron flux calculatedusing 40 precursors clearly diverges from the previous results. The reason for thisbehaviour, is because by removing the 10 most important precursors, the number ofdelayed neutrons emitted decreased, thus the period of the fissile system increasedas explained in Sec. 2.2.1.This deviation from criticality for the case when 40 precursors were used canbe quantified by calculating the fitted reactivity. Indeed, the value of this parameter N o r m a li z ed neu t r on f l u x Figure 4.23: Time evolution for the neutron flux in a water moderated box made ofpure
U, simulated with OpenMC. When the 6 groups are used results shown inblue. Results obtained when 50 individual precursors are used are shown in red, andresults obtained when using 40 precursors are shown in green. In this case it canbe seen that the time evolution of the neutron flux calculated using 40 precursorsclearly diverges from the previous results, showing that the neutron flux grows morerapidly in time.for the 40 precursors calculation was ρ = 0 . ρ fit [pcm] 17(368) 35(347) 111(270)Table 4.19: Results obtained for the reactivity of the water moderated energy depen-dent simulated system in a critical configuration using 6-group, 50 individual and 40individual precursor structures. ystem Configuration Precursor Delayed Simulation Wall-clock OpenMC(TD) Compared Calculated Difference Errorstructure neutron energy time time parameter result with resultMonoenergetic fis-sile Subcritical 1-group Monoenergetic 50 s 12 .
35 h ρ = − ρ = − | ∆ ρ | = 0 pcm δρ = 3 pcm (*)Critical 1-group Monoenergetic 25 s 52 .
97 h ρ =13(70) pcm Point kinetics ρ =10(3) pcm | ∆ ρ | =3 pcm δρ =3 pcm (*) Energy dependentU235
Subcritical 1-group χ ( E ) from 100 µ s 44 .
32 h Λ eff =5 . eff =5 . | ∆Λ eff | =0 .
29 ns δ Λ eff =0 .
56 nsJEFF-3.1.1 β eff =648(38) pcm MCNP β eff =644(6) pcm | ∆ β eff | =4 pcm δβ eff =38 pcm1-group ¯ E g =212 .
31 keV 50 µ s 3 .
51 h Λ eff =5 . eff =5 . | ∆Λ eff | =0 .
29 ns δ Λ eff =0 .
42 nsfrom JEFF-3.1.1 β eff =666(34) pcm MCNP β eff =644(6) pcm | ∆ β eff | =22 pcm δβ eff =34 pcm1-group ¯ E g =501 .
31 keV 50 µ s 3 .
49 h Λ eff =5 . eff =5 . | ∆Λ eff | =0 .
21 ns δ Λ eff =0 .
52 nsfrom ENDF/B-VIII.0 β eff =602(36) pcm MCNP β eff =644(6) pcm | ∆ β eff | =42 pcm δβ eff =36 pcm8-group χ ( E ) from 50 µ s 4 .
32 h Λ eff =5 . eff =5 . | ∆Λ eff | =0 .
29 ns δ Λ eff =0 .
45 nsJEFF-3.1.1 β eff =660(60) pcm MCNP β eff =644(6) pcm | ∆ β eff | =16 pcm δβ eff =60 pcm6-group ¯ E i from 50 µ s 4 .
27 h Λ eff =5 . eff =5 . | ∆Λ eff | =0 .
06 ns δ Λ eff =0 .
29 nsENDF/B-VIII.0 β eff =602(57) pcm MCNP β eff =644(6) pcm | ∆ β eff | =22 pcm δβ eff =57 pcm50 individual ¯ E i from 50 µ s 6 .
43 h Λ eff =5 . eff =5 . | ∆Λ eff | = 0 .
29 ns δ Λ eff =0 .
31 nsENDF/B-VIII.0 β eff =602(57) pcm MCNP β eff =644(6) pcm | ∆ β eff | =42 pcm δβ eff =57 pcmSupercritical 1-group χ ( E ) from 100 µ s 52 .
45 h Λ eff =5 . eff =6 . | ∆Λ eff | =0 .
55 ns δ Λ eff =0 .
29 nsJEFF-3.1.1 β eff =666(56) pcm MCNP β eff =651(6) pcm | ∆ β eff | =15 pcm δβ eff =56 pcm1-group ¯ E g =212 .
31 keV 50 µ s 7 .
84 h Λ eff =5 . eff =6 . | ∆Λ eff | =0 .
55 ns δ Λ eff =0 .
57 nsfrom JEFF-3.1.1 β eff =666(63) pcm MCNP β eff =651(6) pcm | ∆ β eff | =15 pcm δβ eff =63 pcm1-group ¯ E g =501 .
31 keV 50 µ s 7 .
81 h Λ eff =5 . eff =6 . | ∆Λ eff | =0 .
55 ns δ Λ eff =0 .
57 nsfrom ENDF/B-VIII.0 β eff =637(35) pcm MCNP β eff =651(6) pcm | ∆ β eff | =14 pcm δβ eff =35 pcm8-group χ ( E ) from 50 µ s 11 .
19 h Λ eff =6 . eff =6 . | ∆Λ eff | =0 .
03 ns δ Λ eff =0 .
43 nsJEFF-3.1.1 β eff =665(56) pcm MCNP β eff =651(6) pcm | ∆ β eff | =14 pcm δβ eff =56 pcm6-group ¯ E i from 50 µ s 11 .
03 h Λ eff =5 . eff =6 . | ∆Λ eff | =0 .
55 ns δ Λ eff =0 .
57 nsENDF/B-VIII.0 β eff =635(38) pcm MCNP β eff =651(6) pcm | ∆ β eff | =16 pcm δβ eff =38 pcm50 individual ¯ E i from 50 µ s 17 .
24 h Λ eff =5 . eff =6 . | ∆Λ eff | =0 .
55 ns δ Λ eff =0 .
49 nsENDF/B-VIII.0 β eff =621(36) pcm MCNP β eff =651(6) pcm | ∆ β eff | =30 pcm δβ eff =36 pcm Light-water moder-ated energy depen-dent U235
Critical 6-group ¯ E i from 4 s 260 .
05 h ρ =17(368) pcm OpenMC ρ =25(3) pcm | ∆ ρ | =8 pcm δρ =3 pcm (*)50 individual ENDF/B-VIII.0 410 .
76 h ρ =35(347) pcm Criticality | ∆ ρ | =10 pcm40 individual 319 .
65 h ρ =111(270) pcm | ∆ ρ | =86 pcm Table 4.20: Summary of all results obtained with OpenMC(TD).(*) Error was obtained only from the calculatedresult using point kinetics equations or criticality calculation using OpenMC, given that OpenMC(TD) error couldbe improved as explained in Sec. 4.4.2.ith these results, OpenMC(TD) code shows its potential as a Monte Carlotool with the capability to explore how precursor data from nuclear databases im-pacts on results obtained in fissile systems. For instance, since the code can use,in principle, an arbitrary precursor structure, it could be studied how the kineticparameters of a given system responds to changes in the cumulative yield, proba-bility of neutron emission, delayed neutron yield or average delayed neutron energyemitted. In that sense OpenMC(TD) could become a reliable tool to prompt newexperimental data on individual β -delayed neutron emitters.Finally, as it was discussed in Sec. 4.4.2, in order to reduce the associateduncertainties from results obtained, increased simulation times would be required.Regarding this, two possible solutions are proposed: i) use of high computing powerto run the simulations: since the code is already parallelized, it would benefit byhaving a greater number of cores available. Of course, this would require access toinfrastucture, such as supercomputer clusters and ii) implement the variance reduc-tion technique known as “implicit fission” in OpenMC(TD), here, the neutron eitherhas a scattering interation or a fission interaction, and the weight of the neutron ismultiplied by the mean number of fission neutrons produced in the event. By usingthis technique, there is no production of new neutrons during fission, thus reducingthe calculation time; this would require modifications and testings of the code, but itwould be feasible and it could positively impact the current calculation times usingthe same infrastucture used in this thesis.85 hapter 5Summary and conclusions The objective of this work was to explicitly include, in a Monte Carlo simulation, thetime dependence related to the β -delayed neutron emission from individual precur-sors. In order to achieve this, a modified version of OpenMC Monte Carlo simulationcode was developed to include transient capabilities in neutron transport and the op-tion to use individual precursors as β -delayed neutron emitters. This code has beennamed OpenMC(TD) or Time-Dependent OpenMC.OpenMC(TD), in addition to the original OpenMC, includes: i) neutrontime labeling and tracking; ii) monitoring of time dependent parameters in the simu-lation such as neutron flux, reaction rates, neutron current, and total neutron popu-lation; iii) simulation time interval division depending on the detail required for thestudied physics case; iv) a new particle called precursor, which is not transportedand acts as a β -delayed neutron emitter; v) individual precursor properties from nu-clear databases such as precursor cumulative neutron yield, delayed neutron emissionprobability, β -delayed decay constant and average number of delayed neutrons pro-duced per fission ν d ; vi) either precursor N -group grouping capabilities or individualprecursor treatment; vii) forced decay of precursor within each time interval; viii)population control at the end of each time interval using the combing method; andix) a transient source routine to initialize transient simulations.To approach the time modelling of neutron transport and interactions ina experimental nuclear reactor, a fissile system was simulated. OpenMC(TD) was86ested in successively complex systems. Different observables such as reactivity ρ ,effective delayed neutron fraction β eff and effective prompt generation time Λ eff ,obtained with OpenMC(TD) were compared with calculated results, either with ex-act point kinetics solutions (1-group, 6-group, 8-group and 50 individual precursorstructure) or asymptotic decay constant α D (6-group and 50 and 40 individual pre-cursor structure). A summary of the OpenMC(TD) results obtained for the systems,configurations and precursor structures studied in this work is shown in Table 4.20.For the monoenergetic system, using the 1-group precursor structure, differ-ences between OpenMC(TD) and the compared results using point kinetics equationswere within the error of the point kinetics result. Nevertheless, large uncertaintieswere obtained for the reactivity of the subcritical and critical configurations, usingOpenMC(TD).For the light-water moderated energy dependent U system, using the 6-group, 50 and 40 individual precursor structure, differences between OpenMC(TD)and the compared results using criticality calculations with the standard 6-groupprecursor structure, were greater than the error of the criticality calculation. Thesimulation time of 4 s was too short to describe the asymptotical critical behaviourof the system, when the time evolution of the neutron flux increases gradually.For the energy dependent
U system, discrepancies were found in the valueobtained for the effective delayed neutron fraction using JEFF-3.1.1 and ENDF-B/VIII.0 nuclear databases, showing the importance of appropriate nuclear data forindividual precursors. In this case, the simulation time was 100 µ s, with a timeinterval of 100 ns, describing neutron flux prompt changes (prompt drop or promptjump for the subcritical or supercritical configuration, respectively) within the first10 µ s. Both total simulation time and time interval chosen for these cases wereadequate to properly describe the transient behaviour of the neutron flux in thesesystems. Results and its errors can be improved in accuracy by running the simu-lation with longer wall-clock times at CSICCIAN cluster; by applying to comput-ing time outside the institution or by implementing an implicit fission scheme inOpenMC(TD). 87esuming the discussion about the possibility of using OpenMC(TD) tosimulate a full system, such as a nuclear reactor core, according to what has beenlearned and developed in this thesis, this would require i) a complete model of thecore geometry materials, its densities, nuclear cross sections, and to replicate thisprocess with another code, such as MCNP, for its subsequent comparison with respectto k eff ; ii) to read a geometry file at the beginning of each time interval, simulatingin this way the insertion or extraction of the control rods ; and iii) a comparisonwith experimental measurements of reactivity changes.The OpenMC(TD) code, developed in this thesis, shows its potential asa Monte Carlo tool with the capability to explore how precursor data from nucleardatabases impacts on results obtained in fissile systems. In that sense, OpenMC(TD)could become a reliable tool to prompt new experimental data on individual β -delayed neutron emitters. This geometry file will contain the control rods positions at different time. uture work Results obtained with OpenMC(TD) could be compared not only to re-sults obtained from other codes, but also with experimental results from transientmeasurements in nuclear reactors. In the current state of the code, results obtainedwith simulation times of the order of tens of seconds, including individual precursors,would require computational times of the order of months. To decrease this time,reduction variance methods additional to forced decay, combing, and implicit absorp-tion need to be implemented. One of the possible reduction variance methods thatcould be implemented, is implicit fission [43]. The implementation of this methodwould allow to increase the simulation time, as well as decreasing the time intervals,reducing the associated uncertainties of the obtained results from the OpenMC(TD)simulation. Some future problems to study could be reactivity insertion in: i) mod-erated energy dependent system with individual precursors, to quantitatively assessthe relative importance of precursors and thus, prompt experimental measurementsof β -delayed emitter nuclei, and then; ii) reactor fuel and core model, in order toobtain its kinetic parameters and compare with reactivity measurements in a regionof the reactor core, using a reactivimeter.There exists other types of time-dependent problems of special interest inreactor physics, such as burn-up fuel calculations. Although this is a time-dependentcalculation, it requires knowledge of the isotopic abundance of fissile material presentin fuel elements during the fuel period of use, which for an experimental nuclearreactor, is of the order of a few years. Another problem to study could be thecoupling of the time evolution of isotopic abundance obtained using reaction ratescalculated with the Bateman equations, with the OpenMC(TD) code, validated withexperimental measurements during the time when the fuel is used.Nonetheless, the study of the inclusion of time dependence in Monte Carlomethods, would allow to explore other problems where the fuel materials and pre-cursors are not fixed in space, but in movement during the operation time of thenuclear reactor, like fuel in IV-th generation nuclear reactors, such as Molten SaltFast Reactors, where the salt (fuel) moves through the circuit in about 4 s [44] andtransient calculations are needed to take into account this circulation.89 ppendix ADelayed neutron group spectra In this appendix the 8-group β -delayed neutron spectra are shown. × Energy (eV) − × d N / d E ( / e V ) xx Mean = 212300.9 eV
Figure A.1: Group 1, β -delayed neutron energy spectrum from JEFF-3.1.1.90
200 400 600 800 1000 1200 1400 1600 1800 × Energy (eV) − × d N / d E ( / e V ) x x Mean = 616990.3 eV
Figure A.2: Group 2, β -delayed neutron energy spectrum from JEFF-3.1.1. × Energy (eV) − × d N / d E ( / e V ) x x Mean = 288645.9 eV
Figure A.3: Group 3, β -delayed neutron energy spectrum from JEFF-3.1.1.91
500 1000 1500 2000 2500 3000 3500 4000 × Energy (eV) − × d N / d E ( / e V ) Mean = 444633.4 eV x x
Figure A.4: Group 4, β -delayed neutron energy spectrum from JEFF-3.1.1. × Energy (eV) − × d N / d E ( / e V ) Mean = 512005.2 eV x x
Figure A.5: Group 5, β -delayed neutron energy spectrum from JEFF-3.1.1.92 × Energy (eV) − × d N / d E ( / e V ) Mean = 482135.6 eV x x
Figure A.6: Group 6, β -delayed neutron energy spectrum from JEFF-3.1.1. × Energy (eV) − × d N / d E ( / e V ) Mean = 533774.4 eV x x
Figure A.7: Group 7 β -delayed neutron energy spectrum from JEFF-3.1.1.93 × Energy (eV) − × d N / d E ( / e V ) Mean = 586182.1 eV x x Figure A.8: Group 8 β -delayed neutron energy spectrum from JEFF-3.1.1.94 ppendix BSolutions of the Point NeutronKinetics Equations for 1-groupprecursor approximation If the 6-group precursor groups are replaced by 1-group precursor with an effectiveyield fraction β and an effective decay constant given by¯ λ = (cid:88) i β i λ i β , (B.1)then the point kinetics equations become dndt = ρ − β Λ n + λC, (B.2)and dCdt = β Λ n − λC. (B.3)The solutions to Eq. (B.2) and Eq. (B.3) are the time evolution of the neutron andprecursor population, n ( t ) and C ( t ), given by n ( t ) = n (cid:20) ρρ − β exp ( α P t ) − βρ − β exp ( α D t ) (cid:21) , (B.4)and C ( t ) = n (cid:20) ρβ ( ρ − β ) exp ( α P t ) − β Λ λ exp ( α D t ) (cid:21) , (B.5)95here the term α P defined as α P = ρ − β Λ , (B.6)is related to the fast readjustement of the prompt neutron population, which happenson the neutron generation timescale, given a change in the reactivity. On the otherhand, the term α D , defined as α D = λρρ − β , (B.7)corresponds in general to the slower change in the neutron population due to thedelayed source of neutrons, characterized by the precursor decay constant λ .96 ppendix CMonoenergetic fissile system with1-group precursor structure The system described in Sec. 4.1 and in Sec. 4.2 of this work consists of a rectangularbox of 10 cm ×
12 cm ×
20 cm (see Fig. C.1) surrounded by vacuum. All neutronshave the same velocity, making this a mono-energetic problem. Total neutron yieldis fixed, the material cross-sections are constant and there is one precursor group.The system parameters are shown in Table C.1.Figure C.1: Box of 10 cm ×
12 cm ×
20 cm simulated in Sec. 4.2 and Sec. 4.3.97 arameter Value β eff . λ (s − ) 0 . ν . t (cm − ) 1 . f (cm − ) 0 . a (cm − ) 0 . s (cm − ) 0 . v (cm / s) 2 . × Table C.1: Material cross sections and parameters of the monoenergetic system.The mean neutron generation time Λ is given byΛ = (Σ f v ν ) − . (C.1)Here, Σ f is the macroscopic fission cross section, v is neutron speed and ν is theaverage number of neutrons produced per fission. In general, this expression forthe neutron generation depends on the energy, but for a monoenergetic system withconstant cross sections, this value is exact.98 ppendix DSummary of the simulationsperformed in this work Simulations presented in this work were run either in the CSICCIAN cluster fromthe Chilean Nuclear Energy Commission, which is comprised of 32 cores of Intel(R)Xeon(R) CPU E5-2640 v2 @ 2 . . ection 4.2: Monoenergetic fissile system with -group precursor struc-ture All simulations presented in Table D.1 were ran in the CSICCIAN cluster.
Configuration Number Number Number Number Time interval Simulation Wall-clockneutrons precursors batches time intervals length ( ms ) time ( s ) time ( h ) Subcritical 10 ×
60 500 100 50 12 . × . ×
25 5000 10 50 66 . Table D.1: Summary of simulation parameters for monoenergetic fissile system insubcritical, critical, and reactivity insertion configurations, using 1-group precursorstructure.
Section 4.3.1: Energy dependent system - Subcritical configuration
All simulations in Table D.2 were simulated in CSICCIAN cluster, except for studyvi, which was simulated in LIN cluster.
Study Number Number Number Number Time interval Simulation Wall-clockneutrons precursors batches time intervals length (n s ) time ( ms ) time ( h ) i 10 ×
22 1000 100 0 . . × .
05 3 . × .
05 3 . × .
05 4 .
32v 10 × .
05 4 . × .
05 6 . Table D.2: Summary of simulation parameters for energy dependent system in asubcritical configuration, using different precursor structures.100 ection 4.3.2: Energy dependent system - Supercritical configuration
All simulations in Table D.3 were simulated in CSICCIAN cluster, except for thestudy vi, which was simulated in LIN cluster.
Study Number Number Number Number Time interval Simulation Wall-clockneutrons precursors batches time intervals length ( s ) time ( ms ) time ( h ) i 10 ×
10 1000 100 0 . . × .
05 7 . × .
05 7 . × .
05 11 .
19v 10 × .
05 11 . × .
05 17 . Table D.3: Summary of simulation parameters for energy dependent system in asupercritical configuration, using different precursor structures.
Section 4.4: Energy-dependent system with individual precursors andneutron moderator
All simulations in Table D.4 were simulated in LIN cluster.
Precursor Number Number Number Number Time interval Simulation Wall-clockstructure neutrons precursors batches time intervals length ( ms ) time ( s ) time ( h ) × . × . × . Table D.4: Summary of simulation parameters for light-water moderated, energydependent system in a critical configuration, using 6-group, 40 individual, and 50individual precursor structures. 101 ppendix EIndividual precursor data
In this appendix the different precursor structures used in this work are presented.6 -group precursor structureGroup λ ( s − ) β i /β ¯ E (eV) . .
038 4003182 0 . .
213 4665423 0 . .
188 4376344 0 .
311 0 .
407 5524285 1 .
397 0 .
128 5132016 3 .
872 0 .
026 535234Table E.1: 6-group precursor structure used in this work.1020 individual precursor structureNumber Z Symbol A λ ( s − ) I i ¯ E (eV) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . individual precursor structureNumber Z Symbol A λ ( s − ) I i ¯ E (eV) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ibliography [1] D.A. Brown et al. ENDF/B-VIII.0: The 8th Major Release of the NuclearReaction Data Library with CIELO-project Cross Sections, New Standards andThermal Scattering Data. Nuclear Data Sheets , 148:1 – 142, 2018. Special Issueon Nuclear Reaction Data.[2] James J. Duderstadt and Louis J. Hamilton.
Nuclear reactor analysis . WileyNew York, 1976.[3] Bart L. Sjenitzer and J. Eduard Hoogenboom. Dynamic Monte Carlo Methodfor Nuclear Reactor Kinetics Calculations.
Nuclear Science and Engineering ,175(1):94–107, 2013.[4] Antonios G. Mylonakis, M. Varvayanni, D.G.E. Grigoriadis, and N. Catsaros.Developing and investigating a pure Monte-Carlo module for transient neutrontransport analysis.
Annals of Nuclear Energy , 104:103 – 112, 2017.[5] Margaux Faucher.
Coupling between Monte Carlo neutron transport andthermal-hydraulics for the simulation of transients due to reactivity insertions .PhD thesis, Universit´e Paris Saclay, October 2019.[6] G.R. Keepin, T.F. Wimett, and R.K. Zeigler. Delayed neutrons from fissionableisotopes of uranium, plutonium and thorium.
Journal of Nuclear Energy (1954) ,6(1):2 – 21, 1957.[7] Paraskevi Dimitriou et al. A new Reference Database for beta-delayed neutrons.
EPJ Web Conf. , 239:04001, 2020. 1068] International Atomic Energy Agency. IAEA CRP on a Reference Databasefor Beta-Delayed Neutron Emission, 2017. .[9] Paul K. Romano et al. OpenMC: A state-of-the-art Monte Carlo code for re-search and development.
Annals of Nuclear Energy , 82:90 – 97, 2015. Joint In-ternational Conference on Supercomputing in Nuclear Applications and MonteCarlo 2013.[10] Mark Kellett, Olivier Kellett, and Robert Mills.
JEFF Report 20: The JEFF-3.1/-3.1.1 Radioactive Decay Data and Fission Yields Sub-libraries . 01 2009.[11] E E Lewis and W F Miller.
Computational methods of neutron transport . Wiley-Interscience, 1984.[12] G I Bell and S Glasstone.
Nuclear reactor theory . Van Nostrand Reinhold, 1970.[13] M. M. Islam and H. H. Knitter. The Energy Spectrum of Prompt Neutronsfrom the Fission of Uranium-235 by 0.40-MeV Neutrons.
Nuclear Science andEngineering , 50(2):108–114, 1973.[14] M. C. Brady.
Evaluation and application of delayed neutron precursor data .PhD thesis, Los Alamos National Laboratory, United States, 4 1989.[15] G. Rudstam et al. Delayed neutron data for the major actinides.
A report by theWorking Party on International Evaluation Co-operation of the NEA NuclearScience Committee , 01 2002.[16] A Goto, Masaki Fujimaki, Tomoki Fujinawa, N Fukunishi, Hiroo Hasebe, Y Hig-urashi, Kumio Ikegami, E Ikezawa, N Inabe, T Kageyama, O Kamigaito,M Kase, M Kidera, S Kohara, M Kobayashi-Komiyama, M Nagase, K Kumagai,T Maie, T Nakagawa, and Yasushige Yano. Commissioning of RIKEN RI BeamFactory.
Proceedings of the 18th International Conference on Cyclotrons andTheir Applications 2007, CYCLOTRONS 2007 , 01 2007.[17] Daniela Foligno.
New evaluation of delayed-neutron data and associated covari-ances . PhD thesis, Aix Marseille Universit´e, 10 2020.10718] Yican Wu.
Nuclear Data Libraries , pages 181–212. Springer Singapore, Singa-pore, 2019.[19] M Herman and Members of the Cross Sections Evaluation Working Group.ENDF-6 Formats Manual Data Formats and Procedures for the Evaluated Nu-clear Data File ENDF/B-VI and ENDF/B-VII. 6 2009.[20] R. J. Tuttle. Delayed-Neutron Data for Reactor-Physics Analysis.
NuclearScience and Engineering , 56(1):37–71, 1975.[21] A. F . Henry. The Application of Reactor Kinetics to the Analysis of Experi-ments.
Nuclear Science and Engineering , 3(1):52–70, 1958.[22] Sedat Goluoglu and H. L. Dodds. A Time-Dependent, Three-Dimensional Neu-tron Transport Methodology.
Nuclear Science and Engineering , 139(3):248–261,2001.[23] S. Goluoglu et al. A deterministic method for transient, three-dimensional neu-tron transport.
Technical Report Yucca Mountain Project, Las Vegas, Nevada ,(TRN: US0702867), 1 1998.[24] C. L. Bentley, R. DeMeglio, M. E. Dunn, S. Goluoglu, K. Norton, R. Pevey,I. Suslov, and H. Dodds. Development of a hybrid stochastic/deterministicmethod for transient, three dimensional neutron transport.
Transactions of theAmerican Nuclear Society , 66:226–227, 1992.[25] M. Shayesteh and M. Shahriari. Calculation of time-dependent neutronic pa-rameters using Monte Carlo method.
Annals of Nuclear Energy , 36(7):901 –909, 2009.[26] B. Boer, D. Lathouwers, J. L. Kloosterman, T. H. J. J. Van Der Hagen, andG. Strydom. Validation of the DALTON-THERMIX Code System with Tran-sient Analyses of the HTR-10 and Application to the PBMR.
Nuclear Technol-ogy , 170(2):306–321, 2010.[27] R. E. Alcouffe and R. S. Baker. PARTISN: A time-dependent, parallel neutralparticle transport code system.
LA-UR-05-3925 Los Alamos National Labora-tory , 2005. 10828] Frank Graziani.
Computational Methods in Transport: Verification and Valida-tion , volume 62. Springer-Verlag Berlin Heidelberg, 2008.[29] Jeffrey S Rosenthal. Parallel computing and Monte Carlo algorithms.
Far EastJournal of Theoretical Statistics , 4(2):207–236, 2000.[30] J. Spanier and E.M. Gelbard.
Monte Carlo Principles and Neutron TransportProblems . Dover Books on Mathematics Series. Dover Publications, 2008.[31] Sandeep Koranne.
Hierarchical Data Format 5 : HDF5 , pages 191–200. SpringerUS, Boston, MA, 2011.[32] Robert Macfarlane et al. The NJOY Nuclear Data Processing System, Version2016.
Technical Report Los Alamos National Laboratory , (TRN: US1701456), 12017.[33] S.P. Bazzana and J.I. M´arquez Dami´an. IEU-COMP-THERM-014, RA-6 Re-actor: Water Reflected, Water Moderated U(19.77) 3Si2-Al Fuel Plates.
Inter-national Handbook of Evaluated Criticality Safety Benchmark Experiments , 2,2010.[34] S.P. Bazzana. Desarrollo, an´alisis y evaluaci´on de experimentos neutr´onicos enel RA-6. Master’s thesis, Instituto Balseiro, Universidad Nacional de Cuyo,Comisi´on Nacional de Energ´ıa At´omica, Argentina, San Carlos de Bariloche,Argentina, 3 2012. (In spanish).[35] John T. Goorley et al. Initial MCNP6 Release Overview - MCNP6 version 1.0.
LA-UR-13-22934 , 2013.[36] Robin Klein Meulekamp and Steven C. van der Marck. Calculating the EffectiveDelayed Neutron Fraction with Monte Carlo.
Nuclear Science and Engineering ,152(2):142–148, 2006.[37] D. Legrady and J. E. Hoogenboom. Scouting the Feasability of Monte CarloReactor Dynamics Simulations.
International Conference on the Physics ofNuclear Reactors , pages 1–5, 2008. 10938] Thomas Booth. A weight (charge) conserving importance-weighted comb forMonte Carlo.
Conference: American Nuclear Society (ANS) Radiation Protec-tion and Shielding Division topical meeting on advancements and applicationsin radiation protection and shielding , (TRN: 96:002172), 04 1996.[39] Cl´elia de Mulatier, Eric Dumonteil, Alberto Rosso, and Andrea Zoia. Thecritical catastrophe revisited.
Journal of Statistical Mechanics: Theory andExperiment , page P08021, 2015.[40] B. E. Simmons and J. S. King. A Pulsed Neutron Technique for ReactivityDetermination.
Nuclear Science and Engineering , 3(5):595–608, 1958.[41] Masao Yamanaka et al. Effective delayed neutron fraction by Rossi-alphamethod in accelerator-driven system experiments with 100 MeV protons at ky-oto university critical assembly.
Journal of Nuclear Science and Technology ,54(3):293–300, 2017.[42]
International Evaluation of Neutron Cross-Section Standards . InternationalAtomic Energy Agency, Vienna, 2007.[43] Bart L. Sjenitzer and J. Eduard Hoogenboom. Variance reduction for fixed-source Monte Carlo calculations in multiplying systems by improving chain-length statistics.
Annals of Nuclear Energy , 38(10):2195 – 2203, 2011.[44] Mariya Brovchenko et al. Preliminary safety calculations to improve the de-sign of Molten Salt Fast Reactor.