ML18312A310: Difference between revisions

From kanterella
Jump to navigation Jump to search
(Created page by program invented by StriderTol)
 
(Created page by program invented by StriderTol)
Line 19: Line 19:


=Text=
=Text=
{{#Wiki_filter:}}
{{#Wiki_filter:Attachment 5 Computation of Delayed Fission Product Gamma Ray Dose Rates from NCSU PULSTAR Reactor Using a Monte Carlo Number Albedo Approach B.E. Hey ABSTRACT Hey, Brit Elkington.
Computation of Delayed Fisson-Product Gamma-Ray Dose Rates From NCSU PULSTAR Reactor Using a Monte Carlo Number Albedo Approach.(Under the direction of Dr. J. M. Doster.) The dose rate due to decaying fission-product gamma-rays escaping the PULSTAR reactor 10 minutes after a loss of water accident was investigated.
A Monte Carlo simulation of the accident was developed using the departmental VAX-11/730 mini-computer.
Dose rates were calculated at several hundred point detector locations in and around the Burlington facility so that isodose lines (lines of constant dose rate) could be established.
Work preliminary to the simulation included the construction of a table of gamma-ray number albedos for use in sampling and dose estimation.
It was found that many areas would be exposed to gamma radiation fields of significant intensity.
The dose rate computed 20 feet directly over and in line-of-sight of the core was 230 +-10 Rem/hr. on the bay floor was The maximum dose rate computed 175 +-12 mRem/hr. Gamma-ray streaming through bay doors and windows resulted in maximum dose rates of 250 +-13 mRem/hr and 175 +-9 mRem/hr in the control room and loading dock respectively.
The exposure occurring outside the reactor building was due primarily to skyshine and caused dose rates on the order of 4 +-0.3 mRem/hr next to the building.
The dose rate calculated for offices adjacent to the reactor bay were found to be negligible due to attenuation through the bay walls. The number albedos as well as the Monte Carlo accident simulations were benchmarked against available experimental and calculated data found in literature.
Calculated total number albedos agreed with those given by Berger and Raso (9) generally to within 5 per cent. Comparisons between calculated dose rates and those measured at the Livermore Pool-Type Reactor (4) generally show agreement to within 50 per cent and less.
COMPUTATION OF DELAYED FISSION-PRODUCT GAMMA-RAY DOSE RATES FROM NCSU PULSTAR REACTOR USING A MONTE CARLO NUMBER ALBEDO APPROACH by BRIT E. HEY A thesis submitted to the Graduate Faculty of North Carolina State University in partial fulfillment of the requirements for the Degree of Master of Science DEPARTMENT OF NUCLEAR ENGINEERING Raleigh 1 9 8 4 APPROVED BY: Committee Page ii BIOGRAPHY Brit Elkington Hey was born in Lansing, Michigan on November 27, 1957. At the age of five his family moved to Fort Worth, Texas where he lived for two years. Again his family moved to Houston, Texas where he graduated from Scarborough Senior High School in 1976. In August 1977 the author entered North Carolina State University
<NCSU> in Raleigh and received a Bachelor of Science degree in Nuclear Engineering in May 1982. In August 1982 he obtained a fellowship and enrolled in the Nuclear Engineering Graduate School at NCSU where he is currently pursuing a Master of Science degree in Nuclear Engineering.
In January 1980 the education program and Light Company during undergraduate in NCSU. author entered the cooperative was employed by Carolina Power and alternate semesters as an He is now employed by the Gulf States Utilities Company located in Beaumont, Texas. The author Gresham. Mrs. is married to the former Teresa Anne Hey is fully employed in the raising of their two children, Heidi and Houston.
Page iii ACKNOWLEDGEMENT The author wishes to express his deepest gratitude to Dr. Joseph Michael Doster, the chairman of his Advisory Committee.
His advice and council given in the preparation of this thesis proved both sound and invaluable.
The author would also like to express his appreciation to his wife, Teresa, for her patience and understanding and to his little daughter, Heidi, for her sacrifices given unknowingly.
Last but not least the author extends a very dear thanks to his friend, Nita Miller, for her support and encouragement during the course of this study.
Page iv LIST OF TABLES
* LIST OF FIGURES INTRODUCTION .
* LITERATURE REVIEW TABLE OF CONTENTS GAMMA-RAY NUMBER AND DOSE ALBEDOS Introduction
* * * * . * * . * * * * * * *
* Definitions
* * . * . * . . * . . * .
* Generation of Albedos Using Monte Carlo * * . * *
* Comparison of Results to Literature
*..*.*.. Use of Gamma-Ray Albedos in Sampling and Estimation
* . . * . * . . . . . . . . .
A MONTE CARLO MODEL OF DELAYED FISSION PRODUCT GAMMA-RAYS
* * * * * * * * * * * * * * *
* Introduction and General Discussion
* *
* Source Intensity and Energy Distribution
* Importance Sampling From Source Angular and Spatial Distributions
* .
* Geometrical Considerations . * *
* Estimating the Dose * *
* BENCHMARKING
*
* A Comparison to a Simple Analytical Solution
* A Comparison to the Livermore Experiment
* RESULTS AND CONCLUSIONS LIST OF REFERENCES
* . . * . . . Page vi 1 6 16 16 18 23 30 33 42 42 49 54 65 68 72 72 79 86 99 LIST OF TABLES 1. NCSU PULSTAR Reactor data . . 2. Coefficients for least square fits to dose albedo parameters
* . . . . . .
* .
* 3. Total number albedo comparison
: 4. Total dose albedo comparison
: 5. Relative dose rate magnitudes for different times after shutdown . . . . . . . .
: 6. Dose rate comparison to analytical results < mRem/hr) . * . . . . . . . . . . . .
: 7. Dose rate comparison to measured values . Page v Page 3 22 32 32 52 78 82 Page vi LIST OF FIGURES Page 1. Reactor building . . . . . . . . . . . . . . . . . 5 2. Geometry depicting particle reflection from a surf ace . . . . . . . . . . . . . . . .
: 3. Geometry used in albedo computations
: 4. Direction cosines and angles of reflection
: 5. Total fission-product gamma-ray activity versus time after shutdown for reactor operating at 18 25 38 1 watt for 1000 hours . . . . . . . .
53 6. Geometry of core and biological shield 7. Relative axial thermal flux distribution . . . 8. Geometry of simplified problem . . . . . . 9. Lawrence Livermore reactor building . . . . . . 10. Computed dose rates on bay floor <mRem/hr>
: 11. Computed dose rates on bay floor using the C-H formula <mRem/hr> . . . . . . . . . . 12. Computed dose rates on reactor platform <Rem/hr> 13. Computed dose rates in control room <mRem/hr> . 14. Computed dose rates in loading dock <mRem/hr> . 15. Computed dose rates on bay roof <Rem/hr> . . . 16. Computed dose rates outside reactor building ( mRem/hr > * * * * * * * * * * * * * * * * * . . . . . . . . 55 61 74 80 92 93 94 95 96 97 98 INTRODUCTION For a pool-type research reactor a sudden loss of water from the tank is one of the more credible accidents (4). The North Carolina State University
<NCSU) PULSTAR Reactor is susceptable to such an accident.
Table 1 lists its operating characteristics.
The loss of water covering the core could occur since the tank liner has penetrations which represent possible failure modes, permitting the draining of the tank into the reactor room basement.
It is estimated that the time required to drain the tank through tank liner penetrations varies from 2.4 minutes for the 12 X 12 inch beam tube to 12.3 minutes for the 6 inch beam tube facility <20). The power density of the core is low enough to allow air cooling by natural convection after a sudden full power shutdown.
Therefore failure of the fuel pin cladding would not result from loss of tank water (20). The problem arises in the external radiation dose resulting from decaying fission-products within the core. Gamma-rays, which contribute more than 97 per cent to the dose (4), scatter within the core, biological shield, building structure, and air at rates high enough to cause significant dose rates in many areas.
Page 2 The objective of this study is to provide information relevant to several basic questions.
They are: What dose rates are to be expected from this loss of water accident?
What areas must be evacuated?
What routes are the best to take? and; At what time can those evacuated areas be reentered?
This study and its recommendations are to be incorporated into the Safety Analysis Report for the NCSU PULSTAR Reactor and become a permanent part thereof. The data most important to the formulation of safety recommendations are the dose rates found in and around the reactor building.
decreasing due As the dose rates are continuously to radioactive decay of the fission products, all values computed in this study correspond to a reference time of 10 minutes after shutdown.
The areas considered are the bay floor, reactor platform, control room, loading dock, offices, bay roof, and outside the reactor building.
These locations are shown graphically on the isometric figure 1 with the exception of the offices. These are located immediately adjacent to and north of the reactor building.
Page 3 Table 1: NCSU PULSTAR Reactor data <20> Power 1 MW steady state, 2200 MW pulse Maximum Thermal 13 steady state 1.69 x 10/6 Flux 3.72 x 10 pulse Fuel Sintered uo 2 pellets 4 per cent enriched Cladding Zircaloy-2 Moderator-Reflector-Light water Coolant Representative Core 359 Kg uo 2 , 12.6 Kg U-235 Excess Reactivity Temperature Coef. Control Rods Core Dimensions Core Volume Fractions Energy Production Shielding Pool Volume Reactor Bay Dimensions 5 per cent dk/k -8 x 10-5 dk/&deg;F (cold,clean, rods out) Ag-In-Cd (80-15-5>
3 shim safety, 1 pulse 15 7/8 x 15 x 25 inches U02 Gap Water Zr-2 Al 0.3823 0.0155 0.4339 0.0802 0.0881 Reactor on 189 hrs/mo in 1983 44 hrs at 1 MW December 1983 High density and ordinary concrete 14,900 gals height 55 ft width 37 ft length 94 ft Page 4 The most accurate way of obtaining the dose rate information is to actually drain the reactor pool and measure the dose. However this being impractical the best approach is a Monte Carlo simulation of the gamma-ray transport following a loss of water accident.
A computer program, referred to as DOSE here, was developed during the course of this study to perform the simulation.
This approach is the only way to mathematically take into account every significant contribution to the dose. One can model accurately using Monte Carlo all the processes taking place. The drawback of the Monte Carlo method is the large amount of machine time required.
In order to decrease that time and simplify the algorithm a computer program, also utilizing Monte Carlo, was developed to generate a table of differential gamma-ray number albedos. Number rather than dose albedos are used so that contributions from scattered gamma-rays through attenuating material (e.g. concrete walls> can be computed as well. Thus the first part of this study deals with how these albedos are generated and applied. primarily with the Monte Carlo The second part deals simulation of decay fission-product gamma-rays and their contributions to the dose in various locations of the NCSU reactor facility and campus. The time dependence of the dose rate is discussed in the section on source intensity.
Page S Figure 1: Reactor building Page 6 LITERATURE REVIEW The literature reviewed in this study cover several areas due to the many different problems encountered in a simulation of this complexity.
The articles referenced here by no means represent a complete review of either the Monte Carlo method or albedo treatment; only those that have here. a direct bearing on To put this study the present work are reviewed in a better perspective the present day Monte Carlo computer codes available for use in photon or photon-neutron transport calculations are also discussed.
In 1964, an experiment was performed by Knezevich et al. <4> at the Lawrence Radiation Laboratory to study the effects of fission-product decay gamma-rays on the dose rate in and around the reactor building.
The pool of water covering the reactor was drained first to a level of one and then zero foot of water over the top of the core. The significant findings were that undue exposure would not occur unless personnel were situated directly over the reactor core, that delayed neutrons contributed less than 3 per cent of the gamma dose and that the first foot of water over the core decreased the dose rates by 70 percent. Due to the similarities between the results of this experiment and those sought by the present work, the Livermore Page 7 experiment is used as a benchmark to the Monte Carlo program developed here. A comparison between the experimental and calculated results is given in a later chapter. Many authors have written on the method of Monte Carlo for use in particle transport calculations.
An easily understood description of the method is found in Schaeffer's "Reactor Shielding for Nuclear Engineers (2)." Also described are the many types of albedos and their use in dose calculations.
For a more detailed discussion of Monte Carlo the book, "Particle-Transport Simulation With the Monte Carlo Method", written by Carter and Cashell (1) provides an excellent starting point for the reader uninitiated in the aspects of Monte Carlo. There are many articles devoted to the study of specific aspects of Monte Carlo problems.
Shin et al. (16) performed measurements of the neutron and gamma-ray streaming in a cavity-duct system and compared these results to an analysis done by an albedo Monte Carlo method. It was found that for ducts of small radius, the calculations overestimated the gamma-ray streaming due to the fact that the albedo data were generated for plane geometry.
For the case of the PULSTAR reactor the biological shield, which can be considered as a large duct, more than adequately satisfies the large radius criteria.
Page 8 One of the most common and reliable variance reduction techniques used in Monte Carlo is that of splitting the particle when it travels from a region of lesser to greater importance.
Lux developed a relatively straightforward method to obtain approximately optimum regionwise importance (17). It might be obvious that scattering in one region is more important than in another but the relative importances needed to obtain an optimum efficiency are usually not easily found. A relationship central to photon transport is that of the Klein-Nishina differential scattering distribution.
Therefore an accurate yet fast method of sampling from this distribution is extremely important for a Monte Carlo photon transport code processing several thousand histories.
Everett and Cashwell (6) developed formulae that approximate the inverse of the cumulative Klein-Nishina probability distribution.
This allows the direct sampling of the scattered energy with a relative error not exceeding 3.2 per cent over the interval of (0.001,100>
MeV of incident photon energies.
Everett et al. <7> also developed a new Monte Carlo method of sampling the Klein-Nishina probability distribution for the same energy range with a relative error not exceeding 2.2 per cent. Blomquist and Gelbard Cll) compare these methods along with Kahn rejection and others to determine an Page 9 optimum technique.
Due to the great differences in machine capabilities however, no method could be claimed superior over the other. Ideally one would select the best technique for his installation.
For shielding and duct problems where backscattered radiation is the dominant factor contributing to the flux, the concept of albedo is an indispensible tool. Investigators of backscattered gamma-rays are numerous.
Hagqmark et al. <25> measured the angular distribution of the dose rate ratio due to backscattered gamma radiation emerging from concrete, aluminum and steel slabs using Co-60 and Cs-137 sources. Comparisons made with two other experiments showed agreement generally to within 20 per cent. However a comparison made to the Chilton-Huddleston CC-H> formula ClO> fitted to Raso's Monte Carlo calculations (8,9> showed that the experimental values were generally 20 to 35 per cent higher than the calculated results. Wells (26), on the other hand, using Monte Carlo to compute differential dose albedos found that for incident energies greater than 2 MeV, his computed dose albedos were smaller than those reported by Raso and those given by the C-H formula. He goes on to say that his calculated values agree with the experimental data of Hagqmark et al. as well as the calculated values given by either the C-H formula or Raso, an apparent discrepancy.
Page 10 The accuracy of the C-H formula has been the focus of some concern in recent years. Better values for the fitted parameters used by the C-H formula have been obtained C23> as more accurate Monte Carlo calculations of differential dose albedos have become available.
Chilton C24) has devised a modified version of the formula but its greater complexity make it somewhat limiting.
The thermal fission reactor is a common source of gamma radiation encountered by the shielding analyst. Even after the reactor is shut down the decay of fission-products continues to generate significant gamma-ray intensities.
Maienschein (5) performed two sets of experiments to study fission-product gamma-ray emissions which are most important for times after fission greater than one second. In the first set, spectra were measured as a function of time after thermal fission of U-235. The second set of experiments was performed to examine the gamma-ray spectra of several fission-product isotopes in a bare critical assembly.
It was found that after 1000 seconds it was possible to predict the gamma-ray decay rate on the basis of beta-decay information available for the individual fission-products.
The total photon-emission rate and the energy release rate for eight different energy regions as a function of decay time after reactor operation for periods from 1 to 1000 hours are given.
Page 11 A large collection of fission-product data is found in the fission-product decay library of the Evaluated Nuclear Data File <ENDF> which is available from and maintained by the National Nuclear Data Library. The coupled nuclide density rate balance equations utilizing this library can be resolved into a collection of linear chains using the CINDER code and then solved analytically and evaluated numerically.
However George et al. (14) found that results for a unit fission impulse can be fit accurately to a linear combination of exponentials which can subsequently be used to generate the source function.
Attention is now given to a few of the computer codes used for simulating noted by Carter and particle transport.
These codes, as Cashwell Cl), do tend to become obsolete in a few years and it may be that some of those referenced fit that category already. The review does however serve to put into perspective the work cited here. One of the largest and most general of the codes available is MORSE-CG <29). It contains more than 100 subroutines, some of which must be user written. The CG stands for combinatorial geometry which is a general purpose geometry package allowing the description of shapes of physical structures by a combination of several basic geometrical shapes. The code solves neutron, photon and coupled neutron-photon transport problems.
Th.rough the use Page 12 of multigroup cross sections the solutions may be obtained in either the forward or adjoint mode. Time dependence for both shielding and criticality problems is provided.
MORSE also has an albedo option available to the user. The three subroutines
*ALBIN, ALBDO and ALBDOE must be supplied by the user. Also needed are the differential albedos as a function of incident energy, incident polar angle, exit energy, exit polar angle and exit azimuthal angle as well as the total albedos along with three types of cumulative distributions similar to the ones discussed later in this study. Subroutine ALBIN reads the data in while ALBDO is used to sample particle direction and energy. ALBDOE employes a statistical estimation approach to calculate the next-flight contribution to the detector.
Three codes developed by Cashwell et al. at the Los Alamos Scientific Laboratory are MCN, MCG and MCP. MCN (33) is a Monte Carlo neutron transport code. It is referenced here because of the similarities between it and its counterpart photon transport codes. All three handle three dimensional geometry and use point cross section data. The logs of cross sections are stored versus log of energy for a rapid log-log interpolation.
Optional standard variance reduction techniques (e.g. Russian roulette, splitting, and path stretching) are built into the codes. Source information may be inserted in complete generality, although certain standard sources are included.
Page 13 The MCG code (32,35) is suitable for solving a wide variety of gamma transport problems.
The physical processes treated are pair production, compton scatter and photoelectric absorption.
The collision routine assumes photons with energy between 1 keV and 100 MeV. MCP (32,35) is also a photon transport code but has a more sophisticated Monte Carlo collision routine for photons of 1 keV to 15 MeV colliding with atoms of Z = 1, 2, .*. , 94 at rest. It takes into account both incoherent and coherent scattering, fluorescent X-ray emission following photo-electric absorption as well as pair production with local emission of annihilation radiation.
One annihilation photon is treated with a weight of 2W. One also has the option of choosing Thomson scattering (i.e. no energy loss) over Klein-Nishina.
Another code much used is the SORS Monte Carlo photon transport code (34). It treats either 2 or 3 dimensional problems with energies ranging from 1 KeV to 15 MeV. The geometry features are not as general as those of the previous codes. Planes only parallel to the XY, XZ and YZ coordinate planes are allowed. Cross sections for 50 elements are obtained from the LRL photon cross section library tape. Linear interpolation is used for Page 14 intermediate energies.
The code provides for treatment of coherent and incoherent scattering, pair production, photoelectric absorption, and a partial treatment of fluorescence radiation
<i.e. K and L shell only>. Variance reduction techniques include splitting and Russian roulette only. The particle is split into 2 particles when going from a zone of lower to higher importance.
The source package includes a wide variety of geometrical shapes. The source energy can be either monoenergetic or given by a prescribed spectrum.
Robinson (30> developed and applied a coupled Monte Carlo neutron-photon transport code to study the generation and effect of a photon flux in a controlled thermo nuclear reactor <CTR>. The 3 dimensional coupled neutron-photon transport code, POGO, was used to calculate the photon flux incident on the CTR plasma, and the Monte Carlo photon transport code, HUCKFINN, was used to calculate the energy robbed from the plasma via inverse compton events. An inverse compton scatter is one in which the photon gains energy at the expense of the energetic electron.
Although applied to a very specific problem here, it is hoped that this code can be applied to other areas involving neutron-photon coupling.
Examples are: analysis of neutron-photon experiments, nuclear weapons fallout, radiation therapy, and fisson reactor shielding.
Page 15 Another code, developed by Union Carbide and only briefly mentioned here, is PHOTRAN which also is a general purpose photon transport program in complex geometry (31>. Other codes such as the ANDY series of Monte Carlo transport programs <36) are designed primarily for specific time-dependent particle and photon transport applications.
This concludes the brief review of Monte Carlo transport codes. No doubt many have been left out and much more could be said about those mentioned.
From a knowledge of all the various Monte Carlo transport codes available one might ask what need there is to write one himself. The codes are not without their disadvantages.
Although highly flexible and representing the state of the art, the mere size of them can limit their use. There is also a significant amount of work involved in learning how to use the codes and in supplying the necessary user written subroutines.
After all obstacles are overcome the end result is that the user knows how to run a particular code but may not know very much about Monte Carlo or the problem he is attempting to solve. There is then much to be gained by one developing his own Monte Carlo program tailored to his specific application.
The work involved could be comparable to the amount that would be spent in learning to use MORSE, MCG or SORS.
Page 16 GAMMA-RAY NUMBER AND DOSE ALBEDOS Introduction Webster's definition of the word albedo is whiteness or reflective power. In using the albedo concept we are treating the interaction of gamma-rays upon a surf ace as a reflection phenomenon.
Theoretically we know that gamma-rays incident upon a surf ace will scatter throughout the medium: however, if the dimensions of the system are large compared to the mean free path of the gamma-rays scattering in the medium, the albedo treatment is useful. As a preliminary step to the Monte Carlo simulation of the PULSTAR reactor a table of differential gamma-ray number albedos for ordinary concrete is generated using the departmental VAX-11/730 mini-computer.
The reasons for generating this table are as follows. First of all the dimensions of the scattering surfaces making up the reactor shield, bay, and office rooms are very large. Therefore the difference in locations between the incident and exit photons can be neglected.
Secondly, the vast majority of the scattering surf aces are constructed of ordinary concrete.
This material is used in virtually every type of fixed shield and its albedo properties have been thoroughly investigated (2). The lower 12 feet of the biological shield is constructed of high density concrete which is a Page 17 higher Z material.
However Berger and Raso (8,21> show that the total energy albedo for normally incident monoenergetic gamma-rays decreases with increasing atomic number. It is the ref ore assumed that the use of number albedos computed for ordinary concrete will give conservative results. Thirdly, there are tremendous savings in machine time when using albedos over conventional methods of tracking and estimation from interaction point to interaction point. Finally, the differential number albedos contain exit energy information which is lost in the use of dose albedos. Knowing exit energy it is then possible to compute the contribution to the dose from higher order reflections coupled transmission through attenuating material.
with There are however some disadvantages in the use of tabulated differential number albedo data. Besides the obvious increase in computer memory needed to store this data there is also a loss of accuracy in having to interpolate between the values of a table. In this study the advantages are thought to outweigh the disadvantages and a compromise is made between accuracy and the size of the table.
Page 18 Definitions The differential gamma-ray number albedo aCE 0 ,B 0 ,E,e,&#xa2;> as defined here is the differential current out per incident current or a(E 0 ,9 0 ,E,B,&#xa2;> = J(E,e,&#xa2;> ( 1) where E 0 is incident particle energy in MeV, 8 0 is incident polar angle, E is exit particle energy, e is exit polar angle, and is exit azimuthal angle. Figure 2 below illustrates the relationship.
E E 0 SLAB Figure 2: Geometry depicting particle reflection from a surf ace Page 19 As written above this albedo is doubly differential in both exit energy E and solid angle defined by e and &#xa2;;. The total gamma-ray number albedo is given by A!E 0 ,9 0) = !oEMAX!oT!/
2 sin91: a!E 0 ,9 0 ,E,9,&#xa2;> x d &#xa2; d 9 dE < 2 > where EMAX is the maximum exit energy possible given E 0 and 9o. The differential solid angle is sine dB d&#xa2;. If we define Ji= cos e then dJl = -sine de and the equation can also be written as = (EMAX {l!IT a<Eo.J.1 0 ,E,&#xb5;,&#xa2;> d&#xa2; dfl dE }o }o -rr ( 3 ) EMAX is found from conservation of energy and momentum.
The energy after a compton scatter is Eo E = 1 + E 0 Cl -cos8 5 > I 0.511 ( 4) where 8 5 is the scattering angle between the original and scattered line-of-flight.
The cosine of the scattering angle is given by the relationship cos 8 5 = sin 8 0 sin 8 cos tP -cos 9 0 cos 8 ( 5)
Page 20 From figure 2 it can be shown that exit energy E will be a maximum when 8 5 is a minimum or &#xa2; = 0 and 8 = 11/2. Hence cos 8 5 = sin 8 0 and EMAX = ( 6) 1 + E 0 (1 -sin 8 0) I 0.511 For problems where the exit gamma-ray energy does not need to be known exactly or it can be sufficiently approximated by say a single compton scatter then the dose albedo aD<E 0 ,e 0 ,e,&#xa2;> provides the greatest advantage.
The one used here is defined as the differential current out <in dose units) per incident current <in dose units> or ( 7) where K<E> is the flux-to-dose conversion factor. The dose albedo defined above is related to the number albedo defined earlier by the expression Chilton and Huddleston
<lO> developed a semi-empirical formula <the C-H formula> for this differential gamma-ray dose albedo for ordinary concrete using the results of Page 21 Raso's Monte Carlo studies (8). The formula given is C<Eo> K<E 0 ,Bs> x 10 26 + C'<E 0> 1 + cosB 0 sec 9 ( 9) where C<E 0> and C'CE 0) are fitted parameters dependent only upon incident energy E 0 and material properties.
KCE 0 ,B.> is the Klein-Nishina differential energy scattering cross section per electron per steradian.
Note that this is not the same as the Klein-Nishina differential scattering cross section but is given by K(E 0 ,9s> = r/12 CE/Eo>3 CE 0/E + EIE 0 -sin 9s > cm2/e-/ster (10) where r 0 = 2.81794 x 10-13 cm is the classical Thomson radius of the electron.
The Klein-Nishina differential scattering cross section is give by the equation KCE 0 ,9 5) = rf/2 <EIE 0>2 <E 0/E + EIE 0 -sin Bs > cm2/e-/ster
<11> A least squares fit was performed on the parameters CCE> and C'<E> given by the nonweighted model of Chilton and Huddleston to obtain C<E> = b + c E -d E 1*1 ( 12) and Page 22 C' ( E) = e + f E-O. 8 (13) Values for the fitting coefficients b, c, d, e, and f are shown below in table 2. Equation 9 is used successfully in the Monte Carlo simulation as a means of estimating the dose from gamma-ray backscattering off ordinary concrete.
As expected it reduces machine time and does not have the disadvantage of requiring a lot of computer memory. Shown later are comparison runs made using either equation 9 or the number albedo table which result in similar values for the dose. Table 2: Coefficients for least square fits to dose albedo parameters standard linear deviation correlation symbol value (per cent) coefficient b 0.01195 0.01 -----c 0.13994 0.01 0.918 d 0.10190 0.01 0.904 e 0.0070179 0.003 -----f 0.0077319 0.001 0.997 Page 23 Generation of Albedos Using Monte Carlo A table of differential gamma-ray number albedos for ordinary concrete Cl8> is generated using the Monte Carlo method. Values for eleven incident energies and eight angles of incidence are computed.
For each of these eighty-eight incident conditions there are six exit energy intervals, six exit polar angle intervals, and six exit azimuthal angle intervals which total up to approximately twenty-thousand entries that are stored and accessed.
Two-thousand histories are followed for each incident energy and angle. The incident energies E 0 chosen are 0.02, 0.05, 0.075, 0.1, 0.2, 0.5, 1. 0, 1. 25, 2.0, 4.0, and 6.0 MeV. Because the total number albedo changes more rapidly at lower incident energies the intervals are smaller at that end. The eight incident angles are based on the cosine of the angle or &#xb5;0 = cos G 0 and divide evenly the range of 0. 0 to 1.0. Actually the grazing angle is represented by 0.001 and not 0.0. There are 6 exit energy groups for each incident energy E 0 and angle 8 0* They divide evenly the energy range of EMIN to EMAX<E 0 ,e 0> where EMIN = 0.01 MeV is the minimum energy below which histories are terminated and EMAX<E 0 ,e 0> is the maximum exit energy possible as defined earlier. Defining EMAX as such keeps to a minimum Page 24 the number of zeroes stored in the table. The 6 exit polar angle groups are evenly spaced in cose between 0.0 and 1.0. The 6 exit azimuthal angle groups are also evenly spaced between 0.0 and n. Due to the symmetry of the albedo problem it is not angle from -'N to n. necessary to define an exit azimuthal Defining the 9 and intervals as above causes each directional group to subtend an equivalent solid angle. By defining the entering and exiting energy and angular groups, with the exception of incident energy, on evenly spaced intervals one need not search the table to find the appropriate value but can analytically compute the appropriate indices defining the group. Figure 3 illustrates the divisions made in the emergent solid angle as used in ALBEDO. The albedo model itself is simply a homogenized semi-infinite concrete slab 100.0 cm thick. In accordance with photon usual practice, emission are coherent ignored scattering along with and multiple fluorescence radiation and bremsstrahlung which are considered to have little significance
<15).
Page 25 x Figure 3: Geometry used in albedo computations As shown in figure 3 the direction of the incident particle-track vector is parallel to the X-Z plane and goes through the center of the coordinate system having direction cosines 0(, /.1, and 1. Initially Q' = sin9 0 , tf = 0.0, and r = -coseo. The gamma-ray is forced to interact within the slab and is weighted by Page 26 w = 1 -exp< -rr, f > (14) where Of = total cross section of concrete in cm2/g; Pe = pathlength to escape the slab; concrete density in g/cm3. The pathlength to the interaction point is then p = -lnCl -r W>
Cl5) where r is a random variable chosen uniformly over the interval (0.0,l.O>.
using the expression Z = z 0 + 0 p A new z coordinate is then computed (16) where z 0 is initially equal to zero, the surface of the slab. Note that the x and y dependence is not treated since we are simulating this process as a reflection.
Next the type of interaction is chosen. If the energy of the gamma-ray is above 1.022 MeV pair production is possible with the corresponding emission of annihiliation photons. Another random number is chosen and if less than the ratio 17Sl(0'5+cp>
the interaction is a compton scatter, otherwise pair production is assumed. The subscripts s and p stand for scatter and pair production respectively.
Absorption of the gamma-ray is not allowed and the particle is Page 27 weighted by W = ( o; +op> I Oj. In the case of pair production the two 0.511 MeV annihiliation photons are assumed to originate isotropically at the point of positron creation and in opposite directions.
Each is tracked individually.
For a compton scatter a new energy is sampled from the probability Klein-Nishina equation 11. distribution differential function formed from the scattering cross section, The two techniques used in the program ALBEDO for sampling scattered energy are Kahn Rejection Cll) and one developed by Everett, Cashwell, and Turner (7). Kahn Rejection is as surmised a rejection technique and the latter is a new method of sampling the Klein-Nishina probability distribution function.
Both work very well and use approximately the same amount of computing time on the VAX-11/730.
For energies chosen uniformly over the interval of 0.01 to 10.0 MeV the latter method returned energies approximately 1 per cent faster than the former. Once a new energy is sampled it is compared to EMIN. If less than EMIN then the history is terminated; otherwise, the cosine of the scattering angle is found using equation 4. From cos the old direction cosines, and a uniform sampling of the rotational angle about the original line-of-flight, new direction cosines are computed and the scattering process repeats itself. On the average, the photon scatters approximately eight times per history.
Page 28 So far the discussion has described the tracking of gamma-rays in the concrete slab. We now turn to the method of estimation.
At each interaction point, whether it be a compton scatter or pair production, the corresponding photon is forced into each of the thirty-six directional groups and weighted accordingly.
Its energy determines which exit energy group it falls into. The program ALBEDO has the ability to break up one directional group into several subgroups of equal solid angle and then estimate to each one of the subgroups.
The input parameter NEVAL determines the number of subintervals.
For NEVAL=l the subinterval equals the interval.
For NEVAL=2 each interval is divided into two subintervals and each group into four subgroups, etc. For problems where the escape probability or energy of the gamma-ray vary widely within a directional group this gives a more accurate value. When multiplied by the solid angle subtended by the group the differential albedo is then integrated over solid angle within the directional group. In each case the gamma-ray is forced to scatter or emit into the center of the subgroup.
The weight of the gamma-ray entering the directional group is then (17) where W 0 = weight of photon after interaction; W 5 = 1 I 411 for pair production or KCE 0 ,8 5) Io; for compton scatter. It should also be noted that only one of the Page 29 annihilation photons contributes to the albedo directly after emission since the other is always 180 degrees out of phase. To find the pathlength to escape and the differential scattering probability W 5 the cosine of the scattering angle must be obtained.
This is found by use of the equation cos e 5 = ex CX' + /.J/.1 1 + o' < 18 > where ()(, /.J, o and O<', IJ', o' are the direction cosines of the particle before and after scatter respectively.
Page 30 Comparison of Results to Literature A great deal of effort has been devoted to the studying of gamma-rays backscattering from and transmitting through concrete.
Berger <9>, Raso <8,9>, and Wells <26) amassed a considerable amount of data on differential and total albedos using the method of Monte Carlo in the early 60's. There have been numerous attempts to fit this differential data to semi-empirical formulae such as the C-H formula <10> or an exponential curve such as the one developed by Haggmark et al. <28). The results from the Monte Carlo albedo calculations generally agree reasonably well with each other and the experimentally derived data. There are many albedo comparisons found in the literature.
It is not however. the purpose of the present paper to improve or disprove the albedo data present in the literature.
The intent is to generate a table of differential albedos for use in computing the dose-rates resulting from fission-product decay. The present total albedo data is compared to the literature as a means of showing agreement between the albedos calculated here and those of the accepted literature.
Tables 3 and 4 below compare both total number albedos and total dose albedos. With the exception of the C-H formula the results usually agree to within 5 per cent.
Page 31 The albedos calculated for this comparison were generated from the tracking of 1000 gamma-rays incident upon a 20 cm thick concrete slab for each incident energy angle combination.
The mean free path of a 2 MeV gamma-ray in concrete is about 10 cm. The fraction of reflected photons at a depth of 2 mean free paths is almost 100 percent (2). The concrete is taken to be that of ordinary concrete with a density of 2.35 g/cm3. The cutoff energy for history termination is 10 keV. The total albedos were found by numerically integrating the differential albedos over 5 exit polar and 5 exit azimuthal angle intervals.
Table 3: Total number albedo comparison 9 Study 2.0 MeV 1. 0 MeV 0.5 MeV 00 p 0.159 0.215 0.271 BR 0.162 0.221 0.268 R 0.164 0.207 0.275 60&deg; p 0.317 0.351 0.420 BR 0.313 0.390 0.414 R 0.316 0.365 0.419 90&deg; p 0.717 0.733 0.752 BR 0.724 0.744 0.734 R
----------Table 4: Total dose albedo comparison e Study 2.0 MeV 1. 0 MeV oo p 0.017 0.037 BR 0.020 0.040 CH 0.032 0.049 60&deg; p 0.057 0.099 BR 0.055 0.099 CH 0.063 0.094 90&deg; p 0.302 0.347 BR 0.303 0.355 CH 0.322 0.388 P Present paper -1000 histories BR Berger and Raso C9) R Raso (8) 0.5 MeV 0.071 0.074 0.076 0.143 0.146 0.136 0.387 0.395 0.456 Page 32 0.2 MeV 0.258 0.285 0.285 0.407 0.409 0.419 0.703 0.703 -----0.2 MeV 0.139 0.138 0.134 0.211 0.220 0.214 0.463 0.470 0.580 CH Chilton and Huddleston ClO> -numerically integrated Numbers unavailable Page 33 Use of Gamma-Rav Albedos in Sampling and Estimation Before any sampling is done a probability density function Cpdf) must be formed.
is given by a(Eo ,&#xb5;0 ,E,&#xb5;,&#xa2;>
The pdf formed from (19) where is the total number albedo defined by equation 3. For simplicity the parameters E., and &#xb5;:.;will be left off and p(E 0 ,fJc*E,&#xb5;,&#xa2;>
will be referred to as pCE,&#xb5;,&#xa2;>.
The definition of pCE,&#xb5;,&#xa2;> dE dfJ d&#xa2; is then the probability that a particle with incident energy Eu and angle e 0 will be reflected into energy dE about E and into solid angle dfJ d&#xa2; about e and &#xa2;. However pCE,&#xb5;,&#xa2;> is a distribution function in three variables so we define the following pdf s p<E,Jl> = f T1 d&#xa2; -T1 = 2 forrpCE,/J,fP>
dlfi pCE> = lo l p!E,/Jl d/J (20) ( 21>
Page 34 ( 22) We therefore sample for exit energy E' by inverting the equation {E' !1 = j O p CE> dE (23) Once E' is known then exit angle 8' or fl' = cose' is found from inverting t = 2 1:*p!E' .&#xb5;1 pCE') f, fl' = O pCfl/E' > dfl (24) After fl' is found along with E' the exit azimuthal angle &#xa2;' is found by inverting 1&#xa2;* 2 0 p C E , , fl' , &#xa2; > d.f.* r = 3 pCE' ,fl') = 1&#xa2;' 2 O pc&#xa2;1&#xb5;',E'>
dtp <25)
Page 35 where r, , r2 and r3 are all independent uniformly distributed random numbers as defined earlier. The weight assigned to the particle is then ACE 0 ,&#xb5;0> or the total number albedo for a photon having incident energy E 0 and angle defined by &#xb5;o = cos eo. This is the method used in the subroutine ALBSAM to sample for E, &#xb5; and &#xa2;. However for purposes of program efficiency and computer memory the following is done. Letting Ej, &#xb5; k and &#xa2;I be the maximum E, &#xb5; and &#xa2; in the energy or angular group designated by the indices j, k and 1 the following function is formed Since C<Ej,&#xb5;k,&#xa2;1>
is actually a computed table the above is accomplished by j k = z L: L: j'=l k'=l where aj 1 k 1/ = value of diff. j = 1 to J the number of k = 1 to K the number of 1 = 1 to L the number of 1 L: l'=l a '1 L 1/1 I I II J number albedo at Ej, /Jk' exit energy groups; exit polar angle groups; azimuthal angle groups; ( 27) and &#xa2;;I; Page 36 6E = constant exit energy interval; 6f1 = constant exit cose interval; 6&#xa2; = constant exit azimuthal angle interval.
With this notation we can represent the total number albedo ACE 0 ,f1 0> as Therefore Cj, k,I I CJ 1 K,L is a cumulative distribution function (cdf> in three variables.
Thus it can be shown that the right hand sides of equations 23, 24 and 25 can be represented by !o"'i p(El dE = = Cj 1 k,L -Cj-1,k, L Cj)K 1 L -Cj-1) K 1 L Cj 1 1t,I -Cj-l)k 1 1 -Cj,k-1 1 1 + Cj-l,k-1, I Cj,k,L -Cj-1,k,L -Cj,k-1,L + Cj-1 1 k-1 1 L (28) (29) (30) To sample exit energy E from this table of C's we use equation 27 to find a j such that rl Cj,K,L/ cJ,K,L Dj ( 31)
Page 37 Linear interpolation then gives r, -Dj _I E = LlE + E. (32) I D* -Dj-1 I where D 0 is defined to be zero always and E 0 is the lower energy bound <EMIN> for the lowest energy group. Similar interpolations are carried out for &#xb5; and &#xa2; as well using equations 29 and 30. The program DOSE uses direction cosines to define the direction of the particle-track vector. Once the reflected direction given by &#xb5; = cosG and &#xa2; is sampled using ALBSAM it is necessary to translate the incident direction cosines ex, /.? and o into the exit direction cosines ex' , ;:;-' and o: This is accomplished in subroutine TRANS. Shown on figure 4 is the case where a particle is incident upon a surface with the Z-axis directed out from the surface and normal to it. Given are incident direction cosines O<, L? and polar angle of reflection e and azimuthal angle &#xa2;. The emergent direction cosines are then given by ex' = sine cos< r;; + (J > /.j' = sine sin< 9' + (J > o * = cos e (33) (34) (35)
Page 38 z E y Figure 4: Direction cosines and angles of reflection where 'f is given by the relationships ex = cos;; cos/&deg; o = -sin;&deg; to obtain Page 39 (36) where S = 1 for /.] 0 and S = -1 for /.J < 0. Using figure 4 to visualize the opposite case where the Z-axis is directed into the surf ace the emergent direction cosines are found by the similarly derived expressions
()( I = sine cos (,9-' -&#xa2;) (37) tJ* = sine sin<JJ-&#xa2;>
(38) o' = -cose (39) where the angle 9 is again given by equation 36. In the case where the X-axis or Y-axis is normal to the reflecting surf ace the same relationships are equivalent expressions.
used to obtain To estimate the scattered dose at a point detector resulting from a gamma-ray incident upon a reflecting surf ace the photon must be weighted by the differential albedo and in the case of a number albedo an exit energy must be sampled. This is accomplished by the use of function subroutine DNALB.
Page 40 Given the incident direction cosines and those exiting in the direction of the point detector the subroutine UNTRAN, which does the reverse of TRANS, returns &#xb5; and &#xa2;. The cosine of the exiting polar angle,&#xb5;, is found by setting it equal to the absolute value of the emergent direction cosine whose axis is normal to the reflecting surface. The angle;? is obtained by use of equation 36 and the azimuthal angle &#xa2; is found from either equation 33 or 37 depending upon whether the Z-axis is directed inward or outward from the reflecting surface as discussed earlier. Knowing E 0 , &#xb5;0 , &#xb5; and &#xa2; the indices m, i, j and k respectively are computed which symbolize the corresponding energy and angular groups. The indices m and represent E 0 and J.1 0 which as stated earlier are left off for purposes of simplicity.
The differential number albedo in solid angle is then found from = CJ,k 1 I -CJ)k-1)1 -CJ)k)l-1 + CJ)k-1) 1-1 2 &#xa5; fl&#xa2; (40) Given the exiting direction defined by &#xb5;' and &#xa2;' exit energy E' is sampled by inverting the equation 
{E' Jo p(E,fl' ,&#xa2;'> dE' r = P<fl', &#xa2;') {E' = }o pCE/fl' ,&#xa2;' > dE where the right hand side is represented by = Cj,k,I -Cj,k-1,I -Cj,k,l-1 + CjJk-IJl-1 cJ,1<,1 -cJ,k-1,1 -cJ,1<,1-1
+ cJ,k-1,1-1 Page 41 (41) (42) Exit energy E' is then interpolated the same way as described by equations 31 and 32 except that D;, is now equal to equation 42.
Page 42 A MONTE CARLO MODEL OF DELAYED FISSION-PRODUCT GAMMA-RAYS Introduction and General Discussion The computer model DOSE, written in FORTRAN IV, simulates delayed fission-product gamma-rays escaping from the NCSU PULSTAR reactor. Due to the complexities of the problem (e.g. three dimensional geometry and multiple scattering>
Monte Carlo is the method used. Although general purpose Monte Carlo photon transport codes are already available it was decided that a program specifically written for this application would execute much faster, take less memory, and provide a better understanding of the process. Because of the first two advantages the nuclear engineering departmental VAX-11/730 mini-computer was used to perform the computations.
It typically took about 5 hours to obtain a complete dose map of the reactor bay using approximately 200 point detectors.
With 30,000 histories the standard deviations of almost all detectors are under 10 per cent. The program DOSE is general enough to allow the user to obtain dose rate information at as many as 300 different locations at once. Standard deviations for every detector are also computed.
or out of the Detectors can be placed anywhere inside reactor building.
This includes the biological shield, control room, reactor bay, reactor Page 43 platform, loading lock, roof and ground. All that is needed are the detector coordinates with respect to the axes shown in figure 1. All parameters necessary to describe the geometry of the reactor building are specified in subroutine BAY. To change the configuration only BAY has to be altered. The dose rate is a function of the individual detector location.
This requires the estimating of dose rate due to gamma-rays scattering in many different locations which in turn requires time. In order to minimize CPU time there are seven estimating options. These are !BAY, !TRANS, ICSCAT, ISSCAT, IDIREC, !SKY and IALD2. These are set either to 0 or 1, 0 meaning yes and 1 options are defined as follows; meaning no. for IBAY=O The the contributions to the dose from gamma-rays scattering off the reacter bay walls <including control room and loading dock, roof and floor) and the exterior of the biological shield are computed, otherwise they are not, for ITRANS=O the dose transmitted through or backscattered from the bay roof is computed using statistical estimation C2> instead of albedos, for ICSCAT=O the dose from scattered radiation within the core is computed, for ISSCAT=O the dose from radiation scattered within the biological shield is computed, for IDIREC=O the direct contribution to the dose is computed, for ISKY=O the dose due to scattering in the Page 44 atmosphere
<skyshine) is computed in subroutine SKYSHI, and for IALD2=0 the dose due to backscattering off concrete is computed using the C-H formula instead of the number albedo table. The latter is accomplished through the use of function subroutine ALD2. These options allow the computation of dose rates in every location using this single program. They also allow determination of the major sources of radiation since each can be computed individually.
The options apply to all detector locations.
The most important part of the model is the source itself. Small changes in the material composition of the PULSTAR core have a relatively large effect on the resulting dose rate. Table 1 gives the volume fractions of the material making up the core. Gamma-ray attenuation and scattering cross sections are taken from Storm and Israel (3). The core is homogenized and buildup is accounted for by following several scattering events inside the core. Buildup in the core generally accounts for 50 per cent of the dose on the bay floor. For each scattering event the probability of escaping the core is computed and the gamma-ray goes through a series of reflections within the biological shield. This is actually an oversized duct problem for which albedos are idealy suited. Using the albedo table the reflections continue until the gamma-ray escapes the shield or its Page 45 weight or energy become so low that it may be terminated.
Once the escaped photon history is completed, the tracking of the gamma-ray resumes at the point before escape from the core. Using subroutine FSCAT the gamma-ray is forced to scatter again within the core until it too is eventually terminated.
On the average 2 to 3 core scatters are followed per history. The biological shield is thick enough around the base and at the top to allow negligible contribution to the dose from transmission through the shield. It is therefore treated as an opaque structure.
However for detector locations at the top of the biological shield or on the bay roof the dose is due almost entirely to direct, core scattered, and shield scattered radiation escaping from the top opening of the shield. Hence the options IDIREC, ICSCAT and ISSCAT are all set to zero in this case. Eventually a gamma-ray escapes through the top of the biological shield and enters the bay. The most likely point of incidence is with the concrete bay roof. For a thick roof ClO.O cm) and detectors located within the reactor building the escape of gamma-rays through the bay ceiling is not treated. The photon is weighted by the total albedo and forced to reflect off the various concrete structures are met. until criteria for particle history termination Subroutine ALBSAM uses the cumulative Page 46 distribution functions formed from the table of differential number albedos to sample for exit direction and energy. For lower bay detector locations there is negligible contribution from radiation penetrating the shield and therefore only IBAY is set to zero since also has negligible effect. For thin roofs <the Lawrence Livermore benchmark case has a roof of 5/16 inch steel plate) there is a significant reduction in the dose rate due to gamma-rays escaping through the ceiling. The albedo concept is then overly conservative in the case of reflected dose and cannot be used at all for transmitted dose. Statistical estimation is used to estimate from the ceiling instead <IBAY and ITRANS=O>
meaning that the dose is estimated from each scatter within the ceiling and from each reflection within the rest of the bay. Although not allowed to escape the building the gamma-ray can still escape back into the bay as it can contribute further to the dose there. The same procedure is followed for detector locations on the roof where the transmitted dose must be computed.
If the detectors are located in such a way that contributions from direct. core scattered.
and shield scattered radiation are possible then IDIREC. ICSCAT. and ISSCAT must be set to zero as well.
Page 47 For detectors located outside the reactor building at ground level the dose is due almost entirely to backscattering in air. The NCSU reactor bay walls are constructed of concrete 12 inches thick which allows negligible transmission.
The concrete ceiling however is catacombed with the thinnest sections being approximately 3 inches thick. A conservative approximation is made in this case by treating the whole roof as if it were 3 inches thick. To estimate the dose the photon is tracked in the same manner as that used to compute the transmitted dose to detectors located on the roof with the exception that each time the gamma-ray is in the direction of escape through the ceiling the probability of escape is computed, a pathlength is sampled, and the gamma-ray is forced to undergo a series of scatters within an infinite atmosphere until termination criteria are met. For this case ISKY is set to zero which means that only scatters occuring in the atmosphere contribute directly to the dose. Ground effects are neglected.
As expected this case uses the greatest amount of CPU time per history. Once the gamma-ray scattering in the air is terminated the tracking resumes at the point just before its escape from the ceiling. Much has been said about termination criteria.
There are three ways in which the tracking of a photon in a particular region is stopped. The first way is to Page 48 terminate tracking of the photon if its energy falls below the cutoff energy of 20 Kev. When this cutoff energy was lowered to 10 KeV a negligible increase in the dose was observed.
The second way involves Russian roulette <2,19>. Each time the photon weight WC decreases it is compared to its preyious weight or the average weight WAVG before the decrease.
If WC is a certain fraction WMIN below WAVG then Russian roulette is used. A random number r is computed and if it is less than WC I CWAVG x WMIN> the photon weight is increased to WC = WC x WAVG x WMIN, otherwise it is terminated.
Function subroutine
!KEEP accomplishes this task. The optimal WMIN was found by trial and error to be approximately 0.05 for DOSE. The third method of termination is used during dose estimation.
If the maximum dose contribution to all detectors falls below 0.1 per cent of that already accumulated during the history, termination occurs and tracking resumes where it left off.
Page 49 Source Intensity and Energy Distribution The gamma activity is directly proportional to the dose rate. It is determined here by the expression given by Glastone <28) where A;. 0.7 p <t-0.2 _ t-0.2> 0 A = gamma activity in curies; P = average power level in watts; t = number of days after shutdown; tc = number of days after startup. (43) This equation is accurate up to a factor of 2 or less from 10 seconds to several weeks. Due to the fact that the PULSTAR has no fixed operating cycle the average power assumed is 1 MW even though in a typical month such as December 1983 the PULS TAR operated at the peak steady state power of 1 MW for only 44 hours. The time after startup is taken to be infinite.
Thus the reactor is assumed to have reached equillibrium power before shutdown.
These are the most conservative assumptions made in the analysis.
The actual fission-product activity could be much lower. The time after shutdown is taken to be 10 minutes. Substituting these values into equation 43 an activity of 1.9 x 10 6 ci is obtained.
from the loss of water Livermore 2 MW reactor was Page 50 As a comparison the activity experiment conducted at the estimated to be 3.0 x 106 ci after integration of the above equation over its power history. This is 20 per cent below the activity calculated for this reactor using the above assumptions.
The expression used for the approximation of the gamma-ray spectrum is taken from Maienschein (5). It takes the form N<E> = 7.4 exp(-1.10 E> I MeV C44> where N<E> is the yield and E is the gamma-ray energy. The constant corresponds to a total fission-product gamma-ray energy release of 5.9 +-0.7 MeV/fission for times after fission from 1 sec to 108 sec and gamma-ray energies above 2.8 MeV, plus 0.6 +-0.6 MeV/fission for energies below 0.28 MeV, plus 0.3 +-0.2 MeV/fission for times less than 1 sec. The exponent is the same as that used for the prompt fission gamma-rays.
Although roughly approximate we are only concerned with the exponent since the equation is to be formed into a pdf. Also this equation was compared to other equilibrium fission-product gamma-ray spectra <13) which provided energy release rates for seven energy groups along with the effective energy of the group. Converting this energy release data into the yield defined above Page 51 showed that equation 44 predicted a much harder gamma spectrum.
It is therefore thought that its use is at least a conservative approximation to the true fission-product gamma spectrum.
To sample from equation 44 we convert it to the pdf pCE) = 1.1 expC-1.1 E> ( 45) and invert its integral to arrive at the simple expression E = EMIN -ln 51 1.1 (46) The albedo table is good for energies up to 6.0 MeV. Energies sampled above this are set equal to 6.0 MeV. Photons having an energy above 6.0 MeV account for less than 0.2 per cent of the total sampled. The analytical model described later showed that the dose on the bay floor for purely ceiling scattered photons between the energies of 6.0 and 10.0 MeV contributed less than 1 per cent of the dose produced by those emitted between 0.01 and 10.0 MeV. In order to estimate the dose rate at times other than the reference time of 10 minutes after shutdown a log-log plot of data measured by Maienschein is provided in figure 5. The figure gives the fission-product gamma-ray activity versus time after shutdown for a reactor which operated for 1000 hours at 1 watt. We are making the conservative Page 52 assumption that the fission product gamma-ray energy spectrum remains constant with increasing time after shutdown.
From a comparison of Maienschein's gamma-ray energy release rate to gamma-ray activity data, both as a function of time after shutdown, it is seen that the energy spectrum becomes softer with time. Table 5 below uses Maienschein's data to give the relative magnitude of the dose rate for different times after shutdown.
Table 5: Relative dose rate magnitudes for different times after shutdown Time after shutdown 100 sec 10 min 1 hr 1 day 1 wk 1 mo 1 yr Relative Dose Rate 1. 3 1. 0 0.73 0.29 0.14 0.059 0.003 Page 53 -I t-t-a: ;. -1010 -I u w lll ->-t--> -t-u a: 10 9 >-a: 0::: I a: a: L!l 10 8 10 2 20 3 20' io 5 20 6 io' TIME AFTER SHUTDOWN CSECl Figure 5: Total fission-product gamma-ray activity versus time after shutdown for reactor operating at 1 watt for 1000 hours Importance Sampling From Source Angular Distributions The fission-product decay photons Page 54 Spatial are emitted isotropically.
However from figure 6 shown below it is seen th.p.t photons emitted in the solid angle subtended by the biological shield opening have a much better chance of escape than those incident upon the shield wall. It is then advantageous to put more importance on the sampling of photons in that direction.
Discrete importance sampling is used to accomplish this. Given the pdf p(x) defined over the interval [a,bJ we form a modified distribution p'(x) = c pCx> I(x) where c is constant and I(x) is the importance function.
The constant c is found from the property of pdf's defined by the equation l = lbp' <x> dx = c lbp<x> I<x> dx (47) Breaking the distribution p'Cx> into discrete intervals, each having importance In the above can be written as 
-I STD. c 0,\, ,,::_ -.. ' <:;' .... . ' <l ' ...: !::.
* BJ.,,!YYTES CONCRETE CONTROL ROD ACTUATORS TANK Z=O -6 , ' 4 ..> ,'; ' --ct_ Figure 6: Geometry of core and biological shield Page 55 Page 56 1 = c l I/ 1 p<xl dx + ... + In 1 p<xl dx x, 6Xn + .*. + IN l p<xl 6xN dx I N = c L: In l p<x> dx n=l 6Xn N l la Xn pCx> dx f.Xn-1 I = c L: In -a p<x> dx n=l N l Pn -Pn-11 = c L: In (48) n=l where 6xn = interval defining the nth region; In = importance value assigned to the nth region; N = total number of regions. The constant c is then given by 1 c = (49) N L In < Pn -P n-I> n=l To sample from p'<x> an mis found such that Then m-1 f, x r = c L I n C P n -P n -J ) + c Im p C x ) dx n=l x,,,_, (50)
Page 57 =ct:: In <Pn -Pn-/l + c Im <P<x> -Pm-/l (5ll Solving for PCx> we obtain P<x> 1 I m-1 ! = ---5-c} In <Pn -Pn-1> +Pm-I cim fi='1 (52) To obtain x one must invert the function P<x> given by P<xl =
dx ( 53) In this case x =&#xb5;, the cosine of the angle with respect to the Z-axis, a = -1, b = 1, and p<&#xb5;> = 112. The value &#xb5; is then given by p = 2 P<p> -1 (54) To account for the biasing the photon is weighted by w p<&#xb5;> = p' <&#xb5;> 1 = (55) c Im The rotational angle is interval [-fT,iTJ since their importance sampling here. Page 58 sampled uniformly on is little to gain the from This method may be used to importance sample source direction assuming the appropriate values of In are known. The importance of the individual directional regions can not be precisely determined until the problem is solved. The following procedure however will yield adequate estimates of these values. The dose and its absolute standard deviation is computed for one detector at each interval of A test case of a few histories C5000) is run and the standard deviations of the dose originating from source photons emitted in each interval is observed.
Since the variance of the total dose is the sum of the variances of the individual doses it can be determined which interval is causing the greatest statistical error. This interval is also most likely to be one contributing significantly to the dose, hence it is an important region. The importance value In is then increased for that region. This method of trial and error may not be the most eloquent but with a little experimentation it works quite well. The standard deviation can be reduced by as much as 20 per cent for the same number of histories.
The importance values In are read in as part of the input data into the array RIMPFCN>.
Page 59 The decay fission-product gamma-ray source distribution is assumed to be identical to the power distribution of the reactor during steady state operation.
Since the PULS TAR is a heterogeneous reactor the distribution is not entirely smooth. The horizontal distribution generally follows that of a cosine with peaks and dips due to the reflector, water channels and control rods between fuel elements.
Fortunately it was found that the dose rate varied little as the horizontal power distribution was changed. It is therefore adequately approximated by the function p( x) = k cos ( ff x ) 2 <H + h> (56) where p<x> =horizontal power distribution in x <and y); k = constant; H = halfwidth of the core; h = extrapolated distance.
The variable x is defined on the interval [-H,HJ and the extrapolated distance is calculated to be 6.7 cm. The dose rate is much more sensitive to the axial power distribution since photons originating in the top part of the core have a much better chance of escaping the shield. The equation used to describe this distribution is taken Page 60 from Tong and Weisman C27) and takes the form ,(rr k" z) (k" z) p<z> = k --cos 2 L L (57) where pCz> = axial power distribution; k', k" = constants; L = half height of core. This equation is a better approximation to an axial power distribution skewed by rod insertion.
The variable z is defined on the interval [-L,LJ. The constant k" is found from the transcendental equation FCN/Z) = Cl.1585) k" I sink" axial peak power = (58) axial average power where FCN/Z) is the axial nuclear peaking factor and is equal to 1.51 for the PULSTAR (20). This gives a value of 1.227 for k". Shown on figure 7 is a plot of equation 57 superimposed on measured axial thermal flux data taken by a student during a laboratory exercise by activation of a copper wire. The equation predicts the thermal flux in the center quite closely but does not predict the rise in the flux at the ends of the core due to the reflector.
This deviation however represents a relatively small area under the curve and is therefore neglected. 
-Xnld 18IX8 lD CT') o CT') I lD I O' I N -0 I U1 w :r: u z -z 0 -t--U1 0 CL ...J CI: -x CI: Figure 7: Relative axial thermal flux distribution Page 61 Page 62 Source positions sampled from the above distributions without biasing yield dose rates with very high variances and an unreasonable amount of machine time must be spent in reducing them to an acceptable level. This is due to the fact that low energy photons <which make up the majority of those sampled) have an average pathlength of only a few centimeters inside the core. Therefore many photons sampled towards the small probability of towards the surface. center of the core will have a very escape compared to those sampled To reduce the variance we wish to make each history contribute as equal a fraction to the dose as possible.
Using z to represent x, y and z and p(z) to represent p(x), pCy> and pCz) we modify the pdf p(z) representing equations 56 or 57 to obtain the modified pdf p'(z) = p(z) exp<-pc<E>
ID -zl> IN where /lc<E> = attenuation coefficient for the core; D = the z Cx or y> coordinate of the core surface the _photon is exiting; N =normalization constant causing .f p'(z) dz= 1. (59)
Page 63 Equation 59 is basically the original pdf modified by the probability to exit the core unscattered.
This modification is idealy suited to the unscattered gamma-rays emitted from the core but also works very well when the scattered component (buildup>
is included.
The weight to the gamma-ray is then W =WC x p(z) I p'(z) where WC is the photon's current weight. It was found that equation 59 had to be used with caution as its misuse could cause higher variances than that produced by the unmodified pdf. Thus the parameter D is found as follows. In order to determine which side of the core the source photon is exiting, its direction and position must be known. Therefore the source position is sampled from the unmodified pdf first. Once D and the direction are determined the same random number is used to sample from the modified pdf to obtain a coordinate which locates the photon emission closer to the surface it is exiting. Because p'Cz) is rather inefficient to sample from using a rejection technique another method is used. Subtracting the random number r from the cdf created from p'(z) one obtains g(z) = fz p'(z') dz' -r = 0 -L (60)
Page 64 Brent's algorithm is then used to find the zero of g(z). This involves a combination of Newton's method along with bisection.
There is no concern with multiple roots since g(z) has only one zero within the interval [-L,LJ. The same method is used to sample from equation 57. The function subroutine RFISM <Root Finding Importance Sampling Method) performes this task. Normally it takes from 3 to 4 iterations to reduce the relative error to less than 0.01 percent.
Page 65 Geometrical Considerations Th.ere are many structures inside the reactor building that are not modeled in DOSE. Some structures, such as the control rod drive mechanisms located above the core, and the ceiling structure, could have a significant effect on the dose rate. Where approximations are made care is taken to keep the model on the conservative side. The control rod machinery for instance, will tend to attenuate many of the photons streaming out through the top of the biological shield and the structure of the bay ceiling will tend to trap photons directed upward and cause less reflection towards the control room and loading dock. Furthermore windows and wood doors are not modeled and scattering in air while in the reactor building is neglected.
One aspect of the geometry which makes it much simpler to model is the fact that there are no curved surf aces and that almost every planar surf ace is perpendicular to one of the coordinate axes. The two exceptions to the latter are the base of the biological shield and the back of the control room. These are approximated by surfaces perpendicular to the axes. This makes it somewhat easier to compute angles of incidence and reflection.
The direction that the wall is facing is kept track of by the variable !WALL corresponding to necessary part UNTRAN. which takes the 6 sides of the input Page 66 on the values of 1 through 6 of a cube. This is a to subroutines TRANS and All data names of structure locations start with a P for position and end in N, S, E, W, T or B for north, south, east, west, top and bottom. For example PBN stands for position bay north which indicates the position of the north bay wall with respect to the Y-axis. This data is all initialized in a block data subroutine.
Modeling the reflecting surfaces is made simpler by the use of the two subroutines BOX and CUBE. Similar subroutines are used in MORSE-CG.
Given a position and direction within a box or room, BOX returns the pathlength to intersection with the wall, the angle of incidence, and the direction with which the wall is facing CIWALL>. If the initial position is located outside of the box then the returned values pertain to the far side of the box. CUBE on the other hand computes the pathlength to intersection with a solid cube in space, the angle of incidence, the indicator
!WALL, and the indicator IHIT which is 0 if the cube is hit and 1 otherwise.
For both subroutines the caller supplies the 6 wall positions of the box or cube.
Page 67 The biological shield itself is modeled quite accurately by a series of cubes surrounding a box representing the pool. The model is perhaps a little more detailed than necessary but it is done so that the effects of shadowing on the reflected gamma-rays can be observed.
There are several modifications that could be made to DOSE and its subroutines to make it much more general (e.g. inclusion of curved surfaces>, but that is not its purpose. It is sufficiently general to allow easy modifications of the model and yet specific enough such that computing time is not wasted in decision making over various and unnecessary options.
Page 68 Estimating the Dose As discussed in the Introduction General Discussion section there are a number of estimating options that can be used to include the contributions to the dose from different sources or compute them by another method. These options and how they are set for the different cases are explained there as well. Also discussed in the section called, Use of Gamma-Ray Albedos in Sampling and Estimation, are the equations used to form the number albedo distribution functions and their representation by discrete tabulated data. This section details the manner in which the dose rates and their associated statistical variances are obtained.
There are three different methods used here to estimate the dose rate from scattered gamma-rays.
The first method uses the number albedo table. This estimates the dose from a gamma-ray reflecting off a concrete structure and perhaps penetrating some barrier material.
Once the appropriate differential number albedo DNALB and the sampled exit energy E are selected as described by equations 40, 41 and 42 the dose rate due to the reflection is computed by D = s WC DNALB exp(-&#xb5;( E> pc&#xb5; ,(n) K< E> I d2 ( 61>
where S =source intensity (photons/sec>;
WC= photon weight before reflection Ccm-2); DNALB = differential number albedo in solid angle; &#xb5;<E> = attenuation coefficient of barrier material Ccm-1>; pc&#xb5;,&#xa2;> = pathlength through barrier (cm); KCE> = flux to dose conversion factor; Page 69 lld 2 = solid angle subtended by point detector of unit area distance d away. Note that S WC is actually the current due to photons having incident angle e 0 on an elemental area of the reflecting surface. It has the units of photons per crn2 per sec. The approximation of l/d2 for the solid angle subtended by a unit area point detector is valid for large d but can cause problems if the point of reflection or scatter is very close to the detector.
To prevent erroneous results all detectors are placed a few centimeters away from the scattering media. The second method of estimation uses the dose albedo as given by the C-H formula, equation 9. Since an exit energy is not known it is not possible to compute the attenuation coefficient.
Therefore the transmitted dose through a barrier can not be computed.
The dose rate is given by Page 70 (62) Note that S WC KCE 0) is the current incident on the reflecting surface in dose units. The third way of obtaining the dose rate is by a method known as statistical estimation.
This technique is used to estimate from the scattered radiation within a scattering medium. It is given by (63) where pCE 0 ,9 5) is the probability per unit solid angle that a gamma-ray having incident energy E 0 will compton scatter through an angle 8 5* This is given by K<E::i, 8 5 > I oS or the differential over the total Klein-Nishina scattering cross section per electron.
Supposing there are J reflections or scatters in this history then the dose rate computed for this history is X = t Dj < 64) J=l Now supposing there are N particle histories generated the total dose rate is given by u = 1 N Page 71 (65) where U is an estimate of the expected value of the dose rate. An estimate of the variance about the expected value is then given by (66) For each detector a sum of the contributions as well as a sum of the squares of the contributions from each history is stored to compute the variance.
Page 72 BENCHMARKING A Comparison to Simple Analytical Solution As a preliminary benchmark for DOSE a simplified problem is formed, namely to compute the dose rate on the bay floor due to multiple scatters in the bay ceiling only. This excludes all secondary radiation exiting from the biological shield. Since almost all of the gamma-rays escaping from the shield unscattered will exit through the top of the core the horizontal source distribution will have little effect on the dose and is therefore neglected.
Using equation 44 for the energy distribution and equation 57 for the axial source distribution we obtain the differential source term in E and z or photons S<E,z> = S p<E> p(z) (67) sec MeV cm where S is the source intensity in photons per second. If the rectangular shield opening is approximated by a circular opening of equal area then the gamma-rays shining on the ceiling is in the shape of the cone illustrated in figure 8. The differential dose incident upon elemental ceiling area dA is then d(z,E 0 ,8 0) where 1 = distance between source emission point to d.A; K<E> = flux-to-dose conversion factor; &#xb5;c<E> = attenuation coefficient of the core; L = half height of core. Page 73 (68) The height, h, of the ceiling with respect to the center of the core is very large compared to the height, 2L, of the core; therefore, 1 can be approximated by and the variation in the angle of incidence eu can be neglected as the source location varies along the Z-axis. One is then able to analytically integrate the above equation over the height of the core to obtain (69) Multiplying dCE 0 ,8 0) by cose 0 we obtain the differential current in dose units incident upon dA. The differential dose at a point detector due to the reflection off of dA is given by I I I I h I g \ I \ m I \ ; ! I \ I I --4----...J d <POINT DETECTOR CORE \_REFLECTING AREA OF CEILING AXIAL SORCE DISTRIBUTION Figure 8: Geometry of simplified problem Page 74 Page 75 (70) where a 0<E 0 ,e 0 ,e,&#xa2;> is given by the C-H formula and m is the distance from dA to the detector.
All that remains to be done is to integrate d<Ec,6 0 ,e,&#xa2;> over incident energy E 0 and the reflecting area N R2. The differential area dA as seen from figure 8 is given by dA = r d&#xa2; dr. The total dose rate at the detector is then D = (EMAX {d+R 2 }o }d-R x d&#xa2; dr dE <71> where the maximum azimuthal angle is a function of r and is given by -I cos (72) The cosine of the scattering angle is needed for tha C-H formula and is given by equation 5 stated here for convenience cos G 5 = sin 9 0 sin 9 cos i -cos .9 0 cos 9 ( 5) where in this problem cos e = g/m; sinG = rim; sin8 0 = s/l; cosec = h/l; and m = V r 2 + g2; s = Vd 2 + r2 -2 d r cos&#xa2;; 1 =vs2 + h 2. Page 76 The integration in equation 71 is carried out numerically using Simpson's rule. A Monte Carlo simulation of the same problem was made using both the C-H formula and the number albedo table. Below is a comparison table of some typical results as compared to the analytical solutions.
Ideally the Monte Carlo run using the C-H formula should agree with the analytical solution exactly. They are indeed within the standard deviation computed.
The results using the number albedo table are from 10 to 25 percent less than those of the other two methods. From table 3 it is seen that values for the dose albedo computed using the C-H formula are also generally higher than the literature.
The inaccuracies of the C-H formula and proposed modifications to it have been studied in the past (24,26). Deviations in this Page 77 differential albedo with other data can be as great as 50 per cent or more. It is therefore not surprizing to find such variation in the results shown below.
Page 78 Table 6: <mRem/hr)
Dose rate comparison to analytical results Bay Control Reactor Method Floor Room Platform Monte Carlo Number .f\lbedo <1000 hist.) 22.3 +-3.6% 31. 8 +-4.1% 62.1 +-3.6% Monte Carlo Dose Albedo <C-H Formula, 1000 hist. ) 25.5 +-3.1% 41. 7 +-3.4% 77.5 +-3.3% Analytical Solution 25.9 42.7 78.9 Page 79 A Comparison to the Livermore Experiment As a final benchmark for DOSE a series of Monte Carlo runs were made those measured conducted at "Loss-of-Water to compare the calculated dose rates with during the loss of water experiment Livermore.
Experiment This at the report, titled the Livermore Pool-Type Reactor," (4) gave dose rates measured at 10 minutes and 1 hour after shutdown.
Measurements were taken with both 0.0 and 1 foot of water over the core. The dose rates quoted here are the ones taken with 0.0 foot of water over the core. Although this experiment is similar to the type of accident modeled for the PULSTAR there are some major differences.
First of all the Livermore reactor is rated at 2 MW and has a uranium fuel enrichment of 90 per cent. Unlike the PULSTAR this core is constructed of plate type fuel elements of uranium aluminimum alloy. The core and fuel element dimensions obtained from the Livermore reactor design data <12) yield weight percents of the materials comprising the core to be 79 percent aluminum, 18 percent water, and 3 percent uranium. There are also significant differences in the geometry of the biological shield and building.
As shown in figure 9 the reactor building 
.,._ _____ sc' dia. BA 'I ROOF ACTUATORS J 11 5 TEEL (j) 3 PLATFORM __ .... 91_6.._'_.7
.... ,,---__, 3 EI (... dia. -I BIO. 5 HIE L D 26 1 6 II -( *-f_-BA Y FL C OR 4 core -0 . --t> *. 60' 90' Page 80 \:3 I 2C ' (&sect;} -------15C' --------"'L-" Figure 9: Lawrence Livermore reactor building Page 81 is a domed cylindrical structure constructed of steel plate. The biological shield is also cylindrical in shape. In order to model the Livermore reactor facility using DOSE the reactor building is approximated by a rectangular steel structure of the same height and floor area. The cylindrical reactor pool is approximated by a rectangular pool of equal height and volume. Because thermal flux distributions for this core were unavailable the shape is assumed to be the same as that used for the PULSTAR. The fission product activity is taken to be 3.0 x 106 ci which was estimated by those conducting the experiment.
Five runs are made using the above activity as the source intensity.
Both the calculated and measured values are listed in table 7. To identify the representative locations of the computed and measured dose rates table 7 lists numbers under the label "location" which correspond to the circled numbers shown on figure 9. The best agreement is obtained on the reactor balcony which is not in direct line-of-sight with the core. Here the values are within 13 per cent of each other. The dose rates on the roof compare somewhat less favorably.
The maximum computed dose rate of 290 +-12 mRem/hr is more than twice the measured value of 105 mRem/hr. The Page 82 Table 7: Dose rate comparison to measured values Measured Location Dose Rate 1 (Reactor Platform) 195 mRem/hr 2a (Bay Roof) 105 mRem/hr 2b 46.5 mRem/hr 2c 0.58 mRem/hr 2d 0.84 mRem/hr 3 (Line-of-Sight) i 4750 Rem/hr ii 2500 rad/hr 4 (Bay Floor) iii 370 mRem/hr 5a (Outside) 80-60 mRem/hr 5b 55-40 mRem/hr 5c 35-30 mRem/hr 5d < 30 mRem/hr i Estimated value iv v Computed Dose Rate 220 +-6% mRem/hr 290 +-4% mRem/hr 24 +-2% mRem/hr 0.27 +-7% mRem/hr 0.67 +-8% mRem/hr 2925 +-4% Rem/hr 2135 +-4% rad/hr 70 +-6% mRem/hr 720 +-8% mRem/hr 72.3 +-18% mRem/hr 33.2 +-17% mRem/hr 18.7 +-15% mRem/hr 12.8 +-15% mRem/hr ii Measured 1 hour after loss of water iii Maximum dose rate measured on bay floor iv Unconservative case (allows for escape through roof> v Conservative case (no gamma-rays escape bay>
Page 83 taking place is probably due to some added attenuation in unmodeled structures such as control rod actuators, structural beams, etc. The calculated dose rate line-of-sight of the core is 20 feet above and in 38 per cent below the 4750 rem/hr estimated in the experiment.
This apparent underestimation needs to be looked at in light of the fact that the quoted value is estimated and not actually measured.
The dose rate measured at this same point 1 hour after shutdown was 2500 rad/hr. If we use table 5 to account for the decrease in activity one obtains for the computed value 2135 +-85 rad/hr which is now only 15 per cent below the measured dose rate. This indicates that the difference might result more from an error in their estimation rather than an error in the calculated value. The dose rates on the bay floor are computed in two ways. The first allows for the fact that many of the gamma-rays escape through the relatively thin steel roof. The second and more conservative way is to assume that there are more structures in the* bay preventing escape than those modeled. The control rod actuators for instance are located directly above the core and will tend to backscatter many of the gamma-rays escaping the biological shield. Thus a conservative value is obtained by that none of the photons escape the building.
assuming The true Page 84 value should lie somewhere in between as evidenced in table 7. The PULSTAR on the other hand has very thick bay walls Cl ft) and the ceiling is considered to be thick also in the case of detector locations within the reactor building.
The computed dose rates should then be conservative in the case of the PULSTAR. As a final comparison dose rates are computed 3 feet above ground level outside of the reactor building.
The effects taken into account are transmission through the bay walls and skyshine.
The detectors are placed from 60 to 150 feet out from the center of the reactor building.
The computed dose rates decrease from 70 to 10 mRem/hr with standard deviations of approximately 15 per cent. The dose rates measured over the same distances decreased from 80 to 30 mRem/hr. The calculated dose rates would probably be much higher if the shape of the actual domed roof were taken into account as some contribution from transmission through the ceiling would result. It must be kept in mind when comparing these values that several assumptions and approximations have been made in this simulation of the Livermore experiment, some of which are not made in the PULSTAR model. The calculated dose rates however generally agree with the measured values to better than 50 per cent. It is therefore concluded that Page 85 the model is indeed a good one and should give even more accurate results for the PULSTAR for which it was written.
Page 86 RESULTS AND CONCLUSIONS The introduction to this paper pointed out the need to answer some basic questions concerning the dose. They basically ask how much, how long, and where. The analysis is made. keeping the above questions in mind. Detectors are located such that isodose lines, lines of a constant dose rate, may be established.
Approximately one point detector is used per square meter of analyzed surface. All units are in either mRem/hr or Rem/hr as indicated.
A standard deviation is computed for each detector as described by equation 66. Those shown in the following figures however are obtained by averaging the individual deviations of the detectors within the area defined by the isodose lines. The first area of interest is the bay floor. The elevation of all detectors is z = 0 or 0.61 meters above the floor. Figure 10 shows the dose rates computed using 30,000 histories and 5 hours of CPU time. As seen from the figure the hottest spots <150 to 175 +-7% mRem/hr> are located east and south of the reactor. This is due to the off-center location of the core within the biological shield and greater distance between the reactor and south bay wall than between the reactor and north wall. The coolest <least exposed> area is oddly enough adjacent to the biological shield. Due to the fact that 75 per cent of Page 87 the dose is due to shine off the ceiling with the other 25 percent originating from the rest of the bay the biological shield shadows the area around it. Thus the north bay door provides the means of access to the bay floor with the least exposure to gamma radiation.
As* a comparison the same analysis is made using the C-H formula as an estimator instead of the number albedo. The results are similar although about 15 percent higher on the average. The hot spots, as shown on figure 11, are predicted to be in the same locations and the same conclusion as that made above is drawn. One might note that the computed errors are slightly lower in this case since there is not the added statistical variance caused by sampling exit energy. Another point of interest is rates predicted on the bay floor. made concerning dose The average background on the floor while the reactor is operating is measured to be approximately 0.2 mRem/hr. According to figure 10 and table 5 it will take approximately 1 year before the maximum computed dose rate of 175 +-12 mRem/hr is reduced to a value of twice the present background.
The dose rates on the reactor platform shown in figure 12 are computed for detectors located 9 meters above the bay floor. The magnitudes as expected are much higher due Page 88 to the reduced distance between the platform and bay ceiling. The elevation of the detectors are approximately head height if one is standing on the platform.
Because much of the interior surf ace of the biological shield is visible from the platform there is a significant contribution to the dose rate due to gamma-ray reflections within the shield. This accounts for the rapidly increasing dose rates as the detectors are placed closer to the edge of the pool. Other than a location in direct line-of-sight with the core the platform offers the area of greatest exposure to the gamma field and is on the order of several Rem/hr. Figure 13 illustrates the expected dose rates found in the control room. The gamma-rays stream in through the two control room windows and door much as the sun shines in a window. The effect is to cause dose rates on the order of 250 +-13 mRem/hr next to the windows and door decreasing rapidly to a minimum of 5 +-2 mRem/hr in the far corner. The detectors are placed 1 meter above the control room floor. Although very high the dose rates will undoubtedly be reduced drastically by the placement of temporary shields (e.g. lead bricks) over the windows and door.
Page 89 A similar effect is seen in the loading dock area as illustrated in figure 14. The projection is a small landing 12 feet above the bay floor and on the same level as the loading dock. The equipment room is located underneath.
The loading dock is a non-radiation area frequently used as a passageway between buildings located north and south. Dose rates as high as 175 +-9 mRem/hr tapering off to a fraction of a mRem/hr are encountered in this area. As indicated previously attenuation through the door opening into the bay is assumed negligible.
Other areas of interest and ones that may have to be reentered are the offices of the Health Physicist and Nuclear Reactor Operations located adjacent to and north of the reactor building.
Due to the attenuation provided by the solid 1 foot thick concrete bay wall the dose rates predicted are on the order of 10-6 +-10% Rem/hr. This is comparable to the present background and is considered negligible.
Dose rates found on the bay roof range from a few Rem/hr when not in line-of-sight of the core to 25 +-2 Rem/hr when directly overhead 16 meters above its midplane.
Values are shown on figure 15. The dose rate is also computed at a location 7 meters above the midplane of the core and in its direct line-of-sight.
This puts the detector on the same level as the top of the biological shield. A value of 230 +-10 Rem/hr is calculated.
Page 90 The final yet very important area considered is outside the reactor building.
As seen in figure 16 the isodose lines form concentric circles around the reactor building with the reactor core at the center. These dose rates are computed using dry air as the scattering medium. Ground effects are neglected.
The dose is strictly due to skyshine as the bay walls are too thick to allow significant transmission.
As illustrated the dose rate assumes a high value of 4 +-0.3 mRem/hr immediately adjacent to the building to a low value of less than 0.3 +-0.02 mRem/hr for distances greater than 80 meters from the core. Assuming that the natural background outside the reactor building is 200 mRem/yr Cthe average natural background for the United States> the dose rate 10 minutes after the accident would be roughly 200 times the background level at a location next to the building.
In conclusion it can be said that the maximum dose rate occurs on the reactor platform.
Depending upon how close one gets the dose rate can range from a fraction to several hundred rem/hr when located directly over the core. The bay roof is also an area that is subject to a high intensity gamma field, but this is not an area normally occupied by personnel.
From calculations based upon the loss of water experiment conducted at the Livermore Page 91 Pool-Type Reactor the PULSTAR Safety Analysis Report <20) estimates the radiation level outside the control room window to be approximately 100 mR/hr with the pool drained. The corresponding exposure calculated here is 350 +-17 mR/hr or 3 and 1/2 times as high. Both calculations assume an power of 1 MW before shutdown.
It must be understood however, that the values reported here are of necessity somewhat conservative.
As stated earlier the fission-product gamma spectra used in this study is somewhat harder than that measured by Maienschein (5). From a quick glance into the reactor bay one will also notice many other attenuating structures Ceg. stairways, actuators, structural beams) besides the biological shield and bay walls. These are left out in the interest of a conservative analysis.
As far as dose rates outside of the reactor building are concerned it should be noted that personnel located inside the building away from bay openings may be less exposed than those standing immediately outside due to the shielding provided by the building against gamma-rays scattering in air above the reactor building.
NORTH--100 -125 a 8 % 125 -150:J:. 6 % SCALE 0 I 2 3 4 5 meters NORTH BAY DOOR 50-75:t 10% 25-50:t 10% Page 92 Figure 10: Computed dose rates on bay floor CmRem/hr)
Page 93 NORTH 100-125 125-150+/-7
% 150-17S:t.
6 96 75-100 :t. 7 % 50-75 :t. 7 96 25 -50.:t. 8 96 SC AL E ------o 2 3 4 5meters Figure 11: Computed dose rates on bay floor using the C-H formula CmRem/hr)
NORTH ,_-----..._
.,,, ,.,,,,. .,,,, .................... ' / ' , ' / .,...,----
..... ,, ' ,, ,. ... ' I / ' \. .,,. l-10z5% , ' .,,. \. I / I I I I I l \ \ \ ' ' \ ' ' ' .. ' ' .. / ... 0.5: l.t.-5 % ,, --0.3 -0.5 :L 5 % SCA L E I ,, / / I ' I I I I / / / 0 I 2 3 4 Smefers Page 94 10-50.z.496 50-230 .t. 4 % Figure 12: <Rem/hr) Computed dose rates on reactor platform 
....... w .. () 0 51 "O c rt ro p. p. 0 t:ll ro "'! Pl rt ro t:ll ..... ::s n 0 ::s rt "'! 0 ....... "'! 0 0 51 DOOR ---NORTH SCALE 0 I 2 CONTROL ROOM WINDOWS 3 ml'fer s Page 96 t.O ... t.O t.O ... " c c C\l I t.O ...... ---C\l -:x: ..... co Q: -c C) t.O
* l Lw -...J \.) U) Figure 14: Computed dose rates in loading dock <mRem/hr>
Page 97 NORTH ) < I :t 10 % SC ALE 0 2 3 4 5mefESrs Figure 15: Computed dose rates on bay roof <Rem/hr>
Page 98 1-2+/- 7% 0.5-1 t 7% EAST BROUGHTON OR. 10 20 Figure 16: Computed dose rates outside reactor building <mRem/hr)
Page 99 LIST OF REFERENCES
: 1. L. L. Carter and E. D. Cashell, Particle Transport Simulation With the Monte Carlo Method, Technical Information Center, Oak Ridge, Tennessee (1975). 2. N. M. Schaeffer, Reactor Shielding for Nuclear Engineers, U. S. Atomic Energy Commission (1973). 3. E. Storm and H. I. Israel, Photon Cross Sections from 0.001 to 100 MeV for Elements 1 Through 100, University of California, Los Alamos, New Mexico (1967). 4. M. Knezevich, R. L. Kathren, 0. K. Helferich and K. R. Kase, Loss-of-Water Experiment at the Livermore Pool-Type Reactor, Health Physics Journal, 11: pp. 481-487 (1965). 5. F. C. Maienschein, Fission-Product Gamma Rays, Engr. 6. Compendium on Radiation Shielding, 1:76-84, Weslay, New York Cl968}. C. J. Everett and E. the Inverse of Distribution, USAEC Scientific Laboratory D. Cashwell, Approximation for the Klein-Nishina Probability Report LA-4448, Los Alamos (1970). 7. C. J. Everett, E. D. Cashwell and G. D. Turner, On a New Method of Sampling the Klein-Nishina Probability Distributions for All Incident Photon Energies Above 1 KeV, USAEC Report LA-4663, Los Alamos Scientific Laboratory (1971). 8. Dominic J. Reflection Nuc. Sci. Raso, Monte Carlo Calculations on the and Transmission of Scattered Gamma Rays, Eng., 17: pp. 411-418 (1963). 9. Martin J. Berger and Dominic J. Raso, Monte Carlo Calculations of Gamma-Ray Backscattering, Radiat. Res., 12: pp. 20-37 (1960). 10. A. B. Chilton and C. M. Huddleston, A Semiempirical Formula for Differential Dose Albedo for Gamma Rays on Concrete, Nuc. Sci. Eng., 17: pp. 419-424 (1963). 11. R. N. Blomquist and E. M. Gelbard, An Assessment of Existing Klein-Nishina Monte Carlo Sampling Methods, Nuc. Sci. Eng., 83: pp. 380-384 (1983). 12. J. B. Radcliffe, Jr. and E. E. Hill, Lawrence Radiation Laboratory (Livermore), Report UCRL-4919 Rev. (1960).
Page 100 13. Reactor Physics Co:nsta_nts, Second Edition, pp. 629-630 (1963). -. . I ---14. D. George, T. R. England, R. J. Labauve and c. W. Maynard, Photon Sources for Shielding Applications, Trans. Arner. Nuc. Soc., 35: pp. 463-464 (1980}. 15. E. M. Fournie and A. B. *Chilton, _Gamma-Ray Buildup Factors for Concrete Slab Shields Under Slant-Incidence Conditions, Nuc. Sci. Eng:., 76: pp. 66-69 ( 1980). 16. Kazuo Shin, Ryuj i Murakami";
Hiroaki -Taniut:hi, _ 'rumonori Hyodo *:and Yoshiaki'
:. Oka, :;Measurements*-
of -Neutron and Garnma:-Ray Streaming
_ in.i *:.:a Ca vi ty:::Duct x: System. -and ran Analysis.
by* : an Albedo Monte Carlo Method, Nuc. Sci. Eng., 81: pp. 161-171 (1982). 17. Ivan Lux, A Handy Method for Approximate-optimization of Splitting in Monte Carlo, Nuc. Sci. Eng., 83: pp. 198-205 (1983). 18. Ainericari National Standa:rd Guidelines
*on the -Nuclear Analysis and Design of Concrete Radiation Shielding for Nuclear Power' *Plants, American Nuclear Society, ANSI/ANS-6.4 Cl977}. 19. American:
National Standard _for Calculation and _Measurement of Dire<:t and Scatte_red Gamma Radiation from LWR Nuclear Plants, American*
: Nuclear '.-Society., ANSI/ANS-6.6_.l (1979}. NCSV PULSTAR Safety Analysis, C_opy -*Number 80, of *Nuclear-
*:NCSU, -Raleigh, North Carolina*
(1967}. 7:J... M* .
and D._ .J !_ **
._ 0&#xa3;trnnia
_ <1960). , . . -J ... ,,: .. ' -... .' '!' 22. C. M*. ,_ Davisson and L. D. Beach, A Monte Carlo Study bf Gamma.Radiation,-Trans.
Amer.
:*391-(1962->.
: 23. A. B. Chilton, C. M. Davisson, and L. A. Beach, Parameters
_-.for .c
* the C-H .. Albedo Formula for Gamma Rays Reflected from Water-, Concrete, Iron, and_ Lead, Trans. Arner. Nuc. _Soc. , 8: pp. 656-657 C 1965)
* 24. A. B. Chilton, A Modified Formula for. Differential Exposure Albedo for Gamma Rays.Reflected from Concrete, Nuc. Sci. Eng., 27: p. 481 (1967).
Page 101 25. L. G. Haggmark, T. H. Jones, N. E. Scofield and W. J. Gurney, Differential Dose-Rate Measurements of Backscattered Gamma-Rays from Concrete, Aluminum and Steel, Nuc. Sci. Eng.,23: pp. 138-149 (1965). 26. M. B. Wells, Differential Dose Albedos for Calculations of Gamma-Ray Reflection from Concrete, USAEC Report RRA-T56, Radiation Research Associates, Inc. (1964). 27. L. S. Tong and J. Weisman, Thermal Analysis of Pressurized Water Reactors, Amer. Nuc. Soc. Pub., La Grange Park, Illinois, p. 255 (1979). 28. S. Glastone, Principles of Nuclear Reactor Engineering, Van Nostrand, Princeton, p. 120 <1955). 29. Wei-Chien Tung, A Guide to General-Purpose Monte Carlo Gamma-Ray Transport Code With INER-0409/MN-80 (1981). Computer Code MORSE-CG:
Multigroup Neutron and Combinatorial Geometry, 30. P. A. Robinson, Jr., The Development and Application of a Coupled Monte Carlo Neutron-Photon Transport Code, USAEC Report UCRL-51234, University of California Lawrence Livermore Laboratory (1972). 31. C. D. Zerby, J. Agresta, et al., PHOTRAN, A General Purpose Photon Transport Program in Complex Geometry, Technical Report AFWL-TR-65171
<Vols. I-IV>, Union Carbide Corporation, Research Institute (1966-1968).
: 32. E. D. Cashwell, J. R. Neergaard, C. J. Everett, R. G. Schrandt, W. M. Taylor and G. D. Turner, Monte Carlo Photon Codes: MCG and MCP, USAEC Report LA-5157-MS, Los Alamos Scientific Laboratory (1973). 33. E. D. Cashwell, J. R. Neergaard, W. H *. Taylor, G. 34. D. Turner, MCN: A Neutron Monte Carlo Code, USAEC Report LA-4751, Los Alamos Scientific Laboratory (1972). J.W. Kimlinger, F. Plechaty, and SORS Monte Carlo Photon-Transport Code USAEC Report UCRL-50358, University Lawrence Radiation Laboratory (1967). J. R. Terrall, for the CDC 6600, of California
: 35. E .. D. Cashwell, C. J. Everett, and W. M. Taylor, A General Monte Carlo Photon Code, in A Review of the Monte Carlo Method for Radiation Transport Calculations}}

Revision as of 16:27, 19 January 2019

Computation of Delayed Fission Product Gamma Ray Dose Rates from Ncsu Pulstar Reactor Using a Monte Carlo Number Albedo Approach - B. E. Hey
ML18312A310
Person / Time
Site: North Carolina State University
Issue date: 12/31/1984
From: Hey B E
North Carolina State University
To:
Office of Nuclear Reactor Regulation
Shared Package
ML18312A302 List:
References
EPID L-2017-RNW-0026
Download: ML18312A310 (110)
v • d • e


Text

Attachment 5 Computation of Delayed Fission Product Gamma Ray Dose Rates from NCSU PULSTAR Reactor Using a Monte Carlo Number Albedo Approach B.E. Hey ABSTRACT Hey, Brit Elkington.

Computation of Delayed Fisson-Product Gamma-Ray Dose Rates From NCSU PULSTAR Reactor Using a Monte Carlo Number Albedo Approach.(Under the direction of Dr. J. M. Doster.) The dose rate due to decaying fission-product gamma-rays escaping the PULSTAR reactor 10 minutes after a loss of water accident was investigated.

A Monte Carlo simulation of the accident was developed using the departmental VAX-11/730 mini-computer.

Dose rates were calculated at several hundred point detector locations in and around the Burlington facility so that isodose lines (lines of constant dose rate) could be established.

Work preliminary to the simulation included the construction of a table of gamma-ray number albedos for use in sampling and dose estimation.

It was found that many areas would be exposed to gamma radiation fields of significant intensity.

The dose rate computed 20 feet directly over and in line-of-sight of the core was 230 +-10 Rem/hr. on the bay floor was The maximum dose rate computed 175 +-12 mRem/hr. Gamma-ray streaming through bay doors and windows resulted in maximum dose rates of 250 +-13 mRem/hr and 175 +-9 mRem/hr in the control room and loading dock respectively.

The exposure occurring outside the reactor building was due primarily to skyshine and caused dose rates on the order of 4 +-0.3 mRem/hr next to the building.

The dose rate calculated for offices adjacent to the reactor bay were found to be negligible due to attenuation through the bay walls. The number albedos as well as the Monte Carlo accident simulations were benchmarked against available experimental and calculated data found in literature.

Calculated total number albedos agreed with those given by Berger and Raso (9) generally to within 5 per cent. Comparisons between calculated dose rates and those measured at the Livermore Pool-Type Reactor (4) generally show agreement to within 50 per cent and less.

COMPUTATION OF DELAYED FISSION-PRODUCT GAMMA-RAY DOSE RATES FROM NCSU PULSTAR REACTOR USING A MONTE CARLO NUMBER ALBEDO APPROACH by BRIT E. HEY A thesis submitted to the Graduate Faculty of North Carolina State University in partial fulfillment of the requirements for the Degree of Master of Science DEPARTMENT OF NUCLEAR ENGINEERING Raleigh 1 9 8 4 APPROVED BY: Committee Page ii BIOGRAPHY Brit Elkington Hey was born in Lansing, Michigan on November 27, 1957. At the age of five his family moved to Fort Worth, Texas where he lived for two years. Again his family moved to Houston, Texas where he graduated from Scarborough Senior High School in 1976. In August 1977 the author entered North Carolina State University

<NCSU> in Raleigh and received a Bachelor of Science degree in Nuclear Engineering in May 1982. In August 1982 he obtained a fellowship and enrolled in the Nuclear Engineering Graduate School at NCSU where he is currently pursuing a Master of Science degree in Nuclear Engineering.

In January 1980 the education program and Light Company during undergraduate in NCSU. author entered the cooperative was employed by Carolina Power and alternate semesters as an He is now employed by the Gulf States Utilities Company located in Beaumont, Texas. The author Gresham. Mrs. is married to the former Teresa Anne Hey is fully employed in the raising of their two children, Heidi and Houston.

Page iii ACKNOWLEDGEMENT The author wishes to express his deepest gratitude to Dr. Joseph Michael Doster, the chairman of his Advisory Committee.

His advice and council given in the preparation of this thesis proved both sound and invaluable.

The author would also like to express his appreciation to his wife, Teresa, for her patience and understanding and to his little daughter, Heidi, for her sacrifices given unknowingly.

Last but not least the author extends a very dear thanks to his friend, Nita Miller, for her support and encouragement during the course of this study.

Page iv LIST OF TABLES

  • LIST OF FIGURES INTRODUCTION .
  • LITERATURE REVIEW TABLE OF CONTENTS GAMMA-RAY NUMBER AND DOSE ALBEDOS Introduction
  • * * * . * * . * * * * * * *
  • Definitions
  • * . * . * . . * . . * .
  • Generation of Albedos Using Monte Carlo * * . * *
  • Comparison of Results to Literature
  • ..*.*.. Use of Gamma-Ray Albedos in Sampling and Estimation
  • . . * . * . . . . . . . . .

A MONTE CARLO MODEL OF DELAYED FISSION PRODUCT GAMMA-RAYS

  • * * * * * * * * * * * * * *
  • Introduction and General Discussion
  • *
  • Source Intensity and Energy Distribution
  • Importance Sampling From Source Angular and Spatial Distributions
  • .
  • Geometrical Considerations . * *
  • Estimating the Dose * *
  • BENCHMARKING
  • A Comparison to a Simple Analytical Solution
  • A Comparison to the Livermore Experiment
  • RESULTS AND CONCLUSIONS LIST OF REFERENCES
  • . . * . . . Page vi 1 6 16 16 18 23 30 33 42 42 49 54 65 68 72 72 79 86 99 LIST OF TABLES 1. NCSU PULSTAR Reactor data . . 2. Coefficients for least square fits to dose albedo parameters
  • . . . . . .
  • .
  • 3. Total number albedo comparison
4. Total dose albedo comparison
5. Relative dose rate magnitudes for different times after shutdown . . . . . . . .
6. Dose rate comparison to analytical results < mRem/hr) . * . . . . . . . . . . . .
7. Dose rate comparison to measured values . Page v Page 3 22 32 32 52 78 82 Page vi LIST OF FIGURES Page 1. Reactor building . . . . . . . . . . . . . . . . . 5 2. Geometry depicting particle reflection from a surf ace . . . . . . . . . . . . . . . .
3. Geometry used in albedo computations
4. Direction cosines and angles of reflection
5. Total fission-product gamma-ray activity versus time after shutdown for reactor operating at 18 25 38 1 watt for 1000 hours0.0116 days <br />0.278 hours <br />0.00165 weeks <br />3.805e-4 months <br /> . . . . . . . .

53 6. Geometry of core and biological shield 7. Relative axial thermal flux distribution . . . 8. Geometry of simplified problem . . . . . . 9. Lawrence Livermore reactor building . . . . . . 10. Computed dose rates on bay floor <mRem/hr>

11. Computed dose rates on bay floor using the C-H formula <mRem/hr> . . . . . . . . . . 12. Computed dose rates on reactor platform <Rem/hr> 13. Computed dose rates in control room <mRem/hr> . 14. Computed dose rates in loading dock <mRem/hr> . 15. Computed dose rates on bay roof <Rem/hr> . . . 16. Computed dose rates outside reactor building ( mRem/hr > * * * * * * * * * * * * * * * * * . . . . . . . . 55 61 74 80 92 93 94 95 96 97 98 INTRODUCTION For a pool-type research reactor a sudden loss of water from the tank is one of the more credible accidents (4). The North Carolina State University

<NCSU) PULSTAR Reactor is susceptable to such an accident.

Table 1 lists its operating characteristics.

The loss of water covering the core could occur since the tank liner has penetrations which represent possible failure modes, permitting the draining of the tank into the reactor room basement.

It is estimated that the time required to drain the tank through tank liner penetrations varies from 2.4 minutes for the 12 X 12 inch beam tube to 12.3 minutes for the 6 inch beam tube facility <20). The power density of the core is low enough to allow air cooling by natural convection after a sudden full power shutdown.

Therefore failure of the fuel pin cladding would not result from loss of tank water (20). The problem arises in the external radiation dose resulting from decaying fission-products within the core. Gamma-rays, which contribute more than 97 per cent to the dose (4), scatter within the core, biological shield, building structure, and air at rates high enough to cause significant dose rates in many areas.

Page 2 The objective of this study is to provide information relevant to several basic questions.

They are: What dose rates are to be expected from this loss of water accident?

What areas must be evacuated?

What routes are the best to take? and; At what time can those evacuated areas be reentered?

This study and its recommendations are to be incorporated into the Safety Analysis Report for the NCSU PULSTAR Reactor and become a permanent part thereof. The data most important to the formulation of safety recommendations are the dose rates found in and around the reactor building.

decreasing due As the dose rates are continuously to radioactive decay of the fission products, all values computed in this study correspond to a reference time of 10 minutes after shutdown.

The areas considered are the bay floor, reactor platform, control room, loading dock, offices, bay roof, and outside the reactor building.

These locations are shown graphically on the isometric figure 1 with the exception of the offices. These are located immediately adjacent to and north of the reactor building.

Page 3 Table 1: NCSU PULSTAR Reactor data <20> Power 1 MW steady state, 2200 MW pulse Maximum Thermal 13 steady state 1.69 x 10/6 Flux 3.72 x 10 pulse Fuel Sintered uo 2 pellets 4 per cent enriched Cladding Zircaloy-2 Moderator-Reflector-Light water Coolant Representative Core 359 Kg uo 2 , 12.6 Kg U-235 Excess Reactivity Temperature Coef. Control Rods Core Dimensions Core Volume Fractions Energy Production Shielding Pool Volume Reactor Bay Dimensions 5 per cent dk/k -8 x 10-5 dk/°F (cold,clean, rods out) Ag-In-Cd (80-15-5>

3 shim safety, 1 pulse 15 7/8 x 15 x 25 inches U02 Gap Water Zr-2 Al 0.3823 0.0155 0.4339 0.0802 0.0881 Reactor on 189 hrs/mo in 1983 44 hrs at 1 MW December 1983 High density and ordinary concrete 14,900 gals height 55 ft width 37 ft length 94 ft Page 4 The most accurate way of obtaining the dose rate information is to actually drain the reactor pool and measure the dose. However this being impractical the best approach is a Monte Carlo simulation of the gamma-ray transport following a loss of water accident.

A computer program, referred to as DOSE here, was developed during the course of this study to perform the simulation.

This approach is the only way to mathematically take into account every significant contribution to the dose. One can model accurately using Monte Carlo all the processes taking place. The drawback of the Monte Carlo method is the large amount of machine time required.

In order to decrease that time and simplify the algorithm a computer program, also utilizing Monte Carlo, was developed to generate a table of differential gamma-ray number albedos. Number rather than dose albedos are used so that contributions from scattered gamma-rays through attenuating material (e.g. concrete walls> can be computed as well. Thus the first part of this study deals with how these albedos are generated and applied. primarily with the Monte Carlo The second part deals simulation of decay fission-product gamma-rays and their contributions to the dose in various locations of the NCSU reactor facility and campus. The time dependence of the dose rate is discussed in the section on source intensity.

Page S Figure 1: Reactor building Page 6 LITERATURE REVIEW The literature reviewed in this study cover several areas due to the many different problems encountered in a simulation of this complexity.

The articles referenced here by no means represent a complete review of either the Monte Carlo method or albedo treatment; only those that have here. a direct bearing on To put this study the present work are reviewed in a better perspective the present day Monte Carlo computer codes available for use in photon or photon-neutron transport calculations are also discussed.

In 1964, an experiment was performed by Knezevich et al. <4> at the Lawrence Radiation Laboratory to study the effects of fission-product decay gamma-rays on the dose rate in and around the reactor building.

The pool of water covering the reactor was drained first to a level of one and then zero foot of water over the top of the core. The significant findings were that undue exposure would not occur unless personnel were situated directly over the reactor core, that delayed neutrons contributed less than 3 per cent of the gamma dose and that the first foot of water over the core decreased the dose rates by 70 percent. Due to the similarities between the results of this experiment and those sought by the present work, the Livermore Page 7 experiment is used as a benchmark to the Monte Carlo program developed here. A comparison between the experimental and calculated results is given in a later chapter. Many authors have written on the method of Monte Carlo for use in particle transport calculations.

An easily understood description of the method is found in Schaeffer's "Reactor Shielding for Nuclear Engineers (2)." Also described are the many types of albedos and their use in dose calculations.

For a more detailed discussion of Monte Carlo the book, "Particle-Transport Simulation With the Monte Carlo Method", written by Carter and Cashell (1) provides an excellent starting point for the reader uninitiated in the aspects of Monte Carlo. There are many articles devoted to the study of specific aspects of Monte Carlo problems.

Shin et al. (16) performed measurements of the neutron and gamma-ray streaming in a cavity-duct system and compared these results to an analysis done by an albedo Monte Carlo method. It was found that for ducts of small radius, the calculations overestimated the gamma-ray streaming due to the fact that the albedo data were generated for plane geometry.

For the case of the PULSTAR reactor the biological shield, which can be considered as a large duct, more than adequately satisfies the large radius criteria.

Page 8 One of the most common and reliable variance reduction techniques used in Monte Carlo is that of splitting the particle when it travels from a region of lesser to greater importance.

Lux developed a relatively straightforward method to obtain approximately optimum regionwise importance (17). It might be obvious that scattering in one region is more important than in another but the relative importances needed to obtain an optimum efficiency are usually not easily found. A relationship central to photon transport is that of the Klein-Nishina differential scattering distribution.

Therefore an accurate yet fast method of sampling from this distribution is extremely important for a Monte Carlo photon transport code processing several thousand histories.

Everett and Cashwell (6) developed formulae that approximate the inverse of the cumulative Klein-Nishina probability distribution.

This allows the direct sampling of the scattered energy with a relative error not exceeding 3.2 per cent over the interval of (0.001,100>

MeV of incident photon energies.

Everett et al. <7> also developed a new Monte Carlo method of sampling the Klein-Nishina probability distribution for the same energy range with a relative error not exceeding 2.2 per cent. Blomquist and Gelbard Cll) compare these methods along with Kahn rejection and others to determine an Page 9 optimum technique.

Due to the great differences in machine capabilities however, no method could be claimed superior over the other. Ideally one would select the best technique for his installation.

For shielding and duct problems where backscattered radiation is the dominant factor contributing to the flux, the concept of albedo is an indispensible tool. Investigators of backscattered gamma-rays are numerous.

Hagqmark et al. <25> measured the angular distribution of the dose rate ratio due to backscattered gamma radiation emerging from concrete, aluminum and steel slabs using Co-60 and Cs-137 sources. Comparisons made with two other experiments showed agreement generally to within 20 per cent. However a comparison made to the Chilton-Huddleston CC-H> formula ClO> fitted to Raso's Monte Carlo calculations (8,9> showed that the experimental values were generally 20 to 35 per cent higher than the calculated results. Wells (26), on the other hand, using Monte Carlo to compute differential dose albedos found that for incident energies greater than 2 MeV, his computed dose albedos were smaller than those reported by Raso and those given by the C-H formula. He goes on to say that his calculated values agree with the experimental data of Hagqmark et al. as well as the calculated values given by either the C-H formula or Raso, an apparent discrepancy.

Page 10 The accuracy of the C-H formula has been the focus of some concern in recent years. Better values for the fitted parameters used by the C-H formula have been obtained C23> as more accurate Monte Carlo calculations of differential dose albedos have become available.

Chilton C24) has devised a modified version of the formula but its greater complexity make it somewhat limiting.

The thermal fission reactor is a common source of gamma radiation encountered by the shielding analyst. Even after the reactor is shut down the decay of fission-products continues to generate significant gamma-ray intensities.

Maienschein (5) performed two sets of experiments to study fission-product gamma-ray emissions which are most important for times after fission greater than one second. In the first set, spectra were measured as a function of time after thermal fission of U-235. The second set of experiments was performed to examine the gamma-ray spectra of several fission-product isotopes in a bare critical assembly.

It was found that after 1000 seconds it was possible to predict the gamma-ray decay rate on the basis of beta-decay information available for the individual fission-products.

The total photon-emission rate and the energy release rate for eight different energy regions as a function of decay time after reactor operation for periods from 1 to 1000 hours0.0116 days <br />0.278 hours <br />0.00165 weeks <br />3.805e-4 months <br /> are given.

Page 11 A large collection of fission-product data is found in the fission-product decay library of the Evaluated Nuclear Data File <ENDF> which is available from and maintained by the National Nuclear Data Library. The coupled nuclide density rate balance equations utilizing this library can be resolved into a collection of linear chains using the CINDER code and then solved analytically and evaluated numerically.

However George et al. (14) found that results for a unit fission impulse can be fit accurately to a linear combination of exponentials which can subsequently be used to generate the source function.

Attention is now given to a few of the computer codes used for simulating noted by Carter and particle transport.

These codes, as Cashwell Cl), do tend to become obsolete in a few years and it may be that some of those referenced fit that category already. The review does however serve to put into perspective the work cited here. One of the largest and most general of the codes available is MORSE-CG <29). It contains more than 100 subroutines, some of which must be user written. The CG stands for combinatorial geometry which is a general purpose geometry package allowing the description of shapes of physical structures by a combination of several basic geometrical shapes. The code solves neutron, photon and coupled neutron-photon transport problems.

Th.rough the use Page 12 of multigroup cross sections the solutions may be obtained in either the forward or adjoint mode. Time dependence for both shielding and criticality problems is provided.

MORSE also has an albedo option available to the user. The three subroutines

  • ALBIN, ALBDO and ALBDOE must be supplied by the user. Also needed are the differential albedos as a function of incident energy, incident polar angle, exit energy, exit polar angle and exit azimuthal angle as well as the total albedos along with three types of cumulative distributions similar to the ones discussed later in this study. Subroutine ALBIN reads the data in while ALBDO is used to sample particle direction and energy. ALBDOE employes a statistical estimation approach to calculate the next-flight contribution to the detector.

Three codes developed by Cashwell et al. at the Los Alamos Scientific Laboratory are MCN, MCG and MCP. MCN (33) is a Monte Carlo neutron transport code. It is referenced here because of the similarities between it and its counterpart photon transport codes. All three handle three dimensional geometry and use point cross section data. The logs of cross sections are stored versus log of energy for a rapid log-log interpolation.

Optional standard variance reduction techniques (e.g. Russian roulette, splitting, and path stretching) are built into the codes. Source information may be inserted in complete generality, although certain standard sources are included.

Page 13 The MCG code (32,35) is suitable for solving a wide variety of gamma transport problems.

The physical processes treated are pair production, compton scatter and photoelectric absorption.

The collision routine assumes photons with energy between 1 keV and 100 MeV. MCP (32,35) is also a photon transport code but has a more sophisticated Monte Carlo collision routine for photons of 1 keV to 15 MeV colliding with atoms of Z = 1, 2, .*. , 94 at rest. It takes into account both incoherent and coherent scattering, fluorescent X-ray emission following photo-electric absorption as well as pair production with local emission of annihilation radiation.

One annihilation photon is treated with a weight of 2W. One also has the option of choosing Thomson scattering (i.e. no energy loss) over Klein-Nishina.

Another code much used is the SORS Monte Carlo photon transport code (34). It treats either 2 or 3 dimensional problems with energies ranging from 1 KeV to 15 MeV. The geometry features are not as general as those of the previous codes. Planes only parallel to the XY, XZ and YZ coordinate planes are allowed. Cross sections for 50 elements are obtained from the LRL photon cross section library tape. Linear interpolation is used for Page 14 intermediate energies.

The code provides for treatment of coherent and incoherent scattering, pair production, photoelectric absorption, and a partial treatment of fluorescence radiation

<i.e. K and L shell only>. Variance reduction techniques include splitting and Russian roulette only. The particle is split into 2 particles when going from a zone of lower to higher importance.

The source package includes a wide variety of geometrical shapes. The source energy can be either monoenergetic or given by a prescribed spectrum.

Robinson (30> developed and applied a coupled Monte Carlo neutron-photon transport code to study the generation and effect of a photon flux in a controlled thermo nuclear reactor <CTR>. The 3 dimensional coupled neutron-photon transport code, POGO, was used to calculate the photon flux incident on the CTR plasma, and the Monte Carlo photon transport code, HUCKFINN, was used to calculate the energy robbed from the plasma via inverse compton events. An inverse compton scatter is one in which the photon gains energy at the expense of the energetic electron.

Although applied to a very specific problem here, it is hoped that this code can be applied to other areas involving neutron-photon coupling.

Examples are: analysis of neutron-photon experiments, nuclear weapons fallout, radiation therapy, and fisson reactor shielding.

Page 15 Another code, developed by Union Carbide and only briefly mentioned here, is PHOTRAN which also is a general purpose photon transport program in complex geometry (31>. Other codes such as the ANDY series of Monte Carlo transport programs <36) are designed primarily for specific time-dependent particle and photon transport applications.

This concludes the brief review of Monte Carlo transport codes. No doubt many have been left out and much more could be said about those mentioned.

From a knowledge of all the various Monte Carlo transport codes available one might ask what need there is to write one himself. The codes are not without their disadvantages.

Although highly flexible and representing the state of the art, the mere size of them can limit their use. There is also a significant amount of work involved in learning how to use the codes and in supplying the necessary user written subroutines.

After all obstacles are overcome the end result is that the user knows how to run a particular code but may not know very much about Monte Carlo or the problem he is attempting to solve. There is then much to be gained by one developing his own Monte Carlo program tailored to his specific application.

The work involved could be comparable to the amount that would be spent in learning to use MORSE, MCG or SORS.

Page 16 GAMMA-RAY NUMBER AND DOSE ALBEDOS Introduction Webster's definition of the word albedo is whiteness or reflective power. In using the albedo concept we are treating the interaction of gamma-rays upon a surf ace as a reflection phenomenon.

Theoretically we know that gamma-rays incident upon a surf ace will scatter throughout the medium: however, if the dimensions of the system are large compared to the mean free path of the gamma-rays scattering in the medium, the albedo treatment is useful. As a preliminary step to the Monte Carlo simulation of the PULSTAR reactor a table of differential gamma-ray number albedos for ordinary concrete is generated using the departmental VAX-11/730 mini-computer.

The reasons for generating this table are as follows. First of all the dimensions of the scattering surfaces making up the reactor shield, bay, and office rooms are very large. Therefore the difference in locations between the incident and exit photons can be neglected.

Secondly, the vast majority of the scattering surf aces are constructed of ordinary concrete.

This material is used in virtually every type of fixed shield and its albedo properties have been thoroughly investigated (2). The lower 12 feet of the biological shield is constructed of high density concrete which is a Page 17 higher Z material.

However Berger and Raso (8,21> show that the total energy albedo for normally incident monoenergetic gamma-rays decreases with increasing atomic number. It is the ref ore assumed that the use of number albedos computed for ordinary concrete will give conservative results. Thirdly, there are tremendous savings in machine time when using albedos over conventional methods of tracking and estimation from interaction point to interaction point. Finally, the differential number albedos contain exit energy information which is lost in the use of dose albedos. Knowing exit energy it is then possible to compute the contribution to the dose from higher order reflections coupled transmission through attenuating material.

with There are however some disadvantages in the use of tabulated differential number albedo data. Besides the obvious increase in computer memory needed to store this data there is also a loss of accuracy in having to interpolate between the values of a table. In this study the advantages are thought to outweigh the disadvantages and a compromise is made between accuracy and the size of the table.

Page 18 Definitions The differential gamma-ray number albedo aCE 0 ,B 0 ,E,e,¢> as defined here is the differential current out per incident current or a(E 0 ,9 0 ,E,B,¢> = J(E,e,¢> ( 1) where E 0 is incident particle energy in MeV, 8 0 is incident polar angle, E is exit particle energy, e is exit polar angle, and is exit azimuthal angle. Figure 2 below illustrates the relationship.

E E 0 SLAB Figure 2: Geometry depicting particle reflection from a surf ace Page 19 As written above this albedo is doubly differential in both exit energy E and solid angle defined by e and ¢;. The total gamma-ray number albedo is given by A!E 0 ,9 0) = !oEMAX!oT!/

2 sin91: a!E 0 ,9 0 ,E,9,¢> x d ¢ d 9 dE < 2 > where EMAX is the maximum exit energy possible given E 0 and 9o. The differential solid angle is sine dB d¢. If we define Ji= cos e then dJl = -sine de and the equation can also be written as = (EMAX {l!IT a<Eo.J.1 0 ,E,µ,¢> d¢ dfl dE }o }o -rr ( 3 ) EMAX is found from conservation of energy and momentum.

The energy after a compton scatter is Eo E = 1 + E 0 Cl -cos8 5 > I 0.511 ( 4) where 8 5 is the scattering angle between the original and scattered line-of-flight.

The cosine of the scattering angle is given by the relationship cos 8 5 = sin 8 0 sin 8 cos tP -cos 9 0 cos 8 ( 5)

Page 20 From figure 2 it can be shown that exit energy E will be a maximum when 8 5 is a minimum or ¢ = 0 and 8 = 11/2. Hence cos 8 5 = sin 8 0 and EMAX = ( 6) 1 + E 0 (1 -sin 8 0) I 0.511 For problems where the exit gamma-ray energy does not need to be known exactly or it can be sufficiently approximated by say a single compton scatter then the dose albedo aD<E 0 ,e 0 ,e,¢> provides the greatest advantage.

The one used here is defined as the differential current out <in dose units) per incident current <in dose units> or ( 7) where K<E> is the flux-to-dose conversion factor. The dose albedo defined above is related to the number albedo defined earlier by the expression Chilton and Huddleston

<lO> developed a semi-empirical formula <the C-H formula> for this differential gamma-ray dose albedo for ordinary concrete using the results of Page 21 Raso's Monte Carlo studies (8). The formula given is C<Eo> K<E 0 ,Bs> x 10 26 + C'<E 0> 1 + cosB 0 sec 9 ( 9) where C<E 0> and C'CE 0) are fitted parameters dependent only upon incident energy E 0 and material properties.

KCE 0 ,B.> is the Klein-Nishina differential energy scattering cross section per electron per steradian.

Note that this is not the same as the Klein-Nishina differential scattering cross section but is given by K(E 0 ,9s> = r/12 CE/Eo>3 CE 0/E + EIE 0 -sin 9s > cm2/e-/ster (10) where r 0 = 2.81794 x 10-13 cm is the classical Thomson radius of the electron.

The Klein-Nishina differential scattering cross section is give by the equation KCE 0 ,9 5) = rf/2 <EIE 0>2 <E 0/E + EIE 0 -sin Bs > cm2/e-/ster

<11> A least squares fit was performed on the parameters CCE> and C'<E> given by the nonweighted model of Chilton and Huddleston to obtain C<E> = b + c E -d E 1*1 ( 12) and Page 22 C' ( E) = e + f E-O. 8 (13) Values for the fitting coefficients b, c, d, e, and f are shown below in table 2. Equation 9 is used successfully in the Monte Carlo simulation as a means of estimating the dose from gamma-ray backscattering off ordinary concrete.

As expected it reduces machine time and does not have the disadvantage of requiring a lot of computer memory. Shown later are comparison runs made using either equation 9 or the number albedo table which result in similar values for the dose. Table 2: Coefficients for least square fits to dose albedo parameters standard linear deviation correlation symbol value (per cent) coefficient b 0.01195 0.01 -----c 0.13994 0.01 0.918 d 0.10190 0.01 0.904 e 0.0070179 0.003 -----f 0.0077319 0.001 0.997 Page 23 Generation of Albedos Using Monte Carlo A table of differential gamma-ray number albedos for ordinary concrete Cl8> is generated using the Monte Carlo method. Values for eleven incident energies and eight angles of incidence are computed.

For each of these eighty-eight incident conditions there are six exit energy intervals, six exit polar angle intervals, and six exit azimuthal angle intervals which total up to approximately twenty-thousand entries that are stored and accessed.

Two-thousand histories are followed for each incident energy and angle. The incident energies E 0 chosen are 0.02, 0.05, 0.075, 0.1, 0.2, 0.5, 1. 0, 1. 25, 2.0, 4.0, and 6.0 MeV. Because the total number albedo changes more rapidly at lower incident energies the intervals are smaller at that end. The eight incident angles are based on the cosine of the angle or µ0 = cos G 0 and divide evenly the range of 0. 0 to 1.0. Actually the grazing angle is represented by 0.001 and not 0.0. There are 6 exit energy groups for each incident energy E 0 and angle 8 0* They divide evenly the energy range of EMIN to EMAX<E 0 ,e 0> where EMIN = 0.01 MeV is the minimum energy below which histories are terminated and EMAX<E 0 ,e 0> is the maximum exit energy possible as defined earlier. Defining EMAX as such keeps to a minimum Page 24 the number of zeroes stored in the table. The 6 exit polar angle groups are evenly spaced in cose between 0.0 and 1.0. The 6 exit azimuthal angle groups are also evenly spaced between 0.0 and n. Due to the symmetry of the albedo problem it is not angle from -'N to n. necessary to define an exit azimuthal Defining the 9 and intervals as above causes each directional group to subtend an equivalent solid angle. By defining the entering and exiting energy and angular groups, with the exception of incident energy, on evenly spaced intervals one need not search the table to find the appropriate value but can analytically compute the appropriate indices defining the group. Figure 3 illustrates the divisions made in the emergent solid angle as used in ALBEDO. The albedo model itself is simply a homogenized semi-infinite concrete slab 100.0 cm thick. In accordance with photon usual practice, emission are coherent ignored scattering along with and multiple fluorescence radiation and bremsstrahlung which are considered to have little significance

<15).

Page 25 x Figure 3: Geometry used in albedo computations As shown in figure 3 the direction of the incident particle-track vector is parallel to the X-Z plane and goes through the center of the coordinate system having direction cosines 0(, /.1, and 1. Initially Q' = sin9 0 , tf = 0.0, and r = -coseo. The gamma-ray is forced to interact within the slab and is weighted by Page 26 w = 1 -exp< -rr, f > (14) where Of = total cross section of concrete in cm2/g; Pe = pathlength to escape the slab; concrete density in g/cm3. The pathlength to the interaction point is then p = -lnCl -r W>

Cl5) where r is a random variable chosen uniformly over the interval (0.0,l.O>.

using the expression Z = z 0 + 0 p A new z coordinate is then computed (16) where z 0 is initially equal to zero, the surface of the slab. Note that the x and y dependence is not treated since we are simulating this process as a reflection.

Next the type of interaction is chosen. If the energy of the gamma-ray is above 1.022 MeV pair production is possible with the corresponding emission of annihiliation photons. Another random number is chosen and if less than the ratio 17Sl(0'5+cp>

the interaction is a compton scatter, otherwise pair production is assumed. The subscripts s and p stand for scatter and pair production respectively.

Absorption of the gamma-ray is not allowed and the particle is Page 27 weighted by W = ( o; +op> I Oj. In the case of pair production the two 0.511 MeV annihiliation photons are assumed to originate isotropically at the point of positron creation and in opposite directions.

Each is tracked individually.

For a compton scatter a new energy is sampled from the probability Klein-Nishina equation 11. distribution differential function formed from the scattering cross section, The two techniques used in the program ALBEDO for sampling scattered energy are Kahn Rejection Cll) and one developed by Everett, Cashwell, and Turner (7). Kahn Rejection is as surmised a rejection technique and the latter is a new method of sampling the Klein-Nishina probability distribution function.

Both work very well and use approximately the same amount of computing time on the VAX-11/730.

For energies chosen uniformly over the interval of 0.01 to 10.0 MeV the latter method returned energies approximately 1 per cent faster than the former. Once a new energy is sampled it is compared to EMIN. If less than EMIN then the history is terminated; otherwise, the cosine of the scattering angle is found using equation 4. From cos the old direction cosines, and a uniform sampling of the rotational angle about the original line-of-flight, new direction cosines are computed and the scattering process repeats itself. On the average, the photon scatters approximately eight times per history.

Page 28 So far the discussion has described the tracking of gamma-rays in the concrete slab. We now turn to the method of estimation.

At each interaction point, whether it be a compton scatter or pair production, the corresponding photon is forced into each of the thirty-six directional groups and weighted accordingly.

Its energy determines which exit energy group it falls into. The program ALBEDO has the ability to break up one directional group into several subgroups of equal solid angle and then estimate to each one of the subgroups.

The input parameter NEVAL determines the number of subintervals.

For NEVAL=l the subinterval equals the interval.

For NEVAL=2 each interval is divided into two subintervals and each group into four subgroups, etc. For problems where the escape probability or energy of the gamma-ray vary widely within a directional group this gives a more accurate value. When multiplied by the solid angle subtended by the group the differential albedo is then integrated over solid angle within the directional group. In each case the gamma-ray is forced to scatter or emit into the center of the subgroup.

The weight of the gamma-ray entering the directional group is then (17) where W 0 = weight of photon after interaction; W 5 = 1 I 411 for pair production or KCE 0 ,8 5) Io; for compton scatter. It should also be noted that only one of the Page 29 annihilation photons contributes to the albedo directly after emission since the other is always 180 degrees out of phase. To find the pathlength to escape and the differential scattering probability W 5 the cosine of the scattering angle must be obtained.

This is found by use of the equation cos e 5 = ex CX' + /.J/.1 1 + o' < 18 > where ()(, /.J, o and O<', IJ', o' are the direction cosines of the particle before and after scatter respectively.

Page 30 Comparison of Results to Literature A great deal of effort has been devoted to the studying of gamma-rays backscattering from and transmitting through concrete.

Berger <9>, Raso <8,9>, and Wells <26) amassed a considerable amount of data on differential and total albedos using the method of Monte Carlo in the early 60's. There have been numerous attempts to fit this differential data to semi-empirical formulae such as the C-H formula <10> or an exponential curve such as the one developed by Haggmark et al. <28). The results from the Monte Carlo albedo calculations generally agree reasonably well with each other and the experimentally derived data. There are many albedo comparisons found in the literature.

It is not however. the purpose of the present paper to improve or disprove the albedo data present in the literature.

The intent is to generate a table of differential albedos for use in computing the dose-rates resulting from fission-product decay. The present total albedo data is compared to the literature as a means of showing agreement between the albedos calculated here and those of the accepted literature.

Tables 3 and 4 below compare both total number albedos and total dose albedos. With the exception of the C-H formula the results usually agree to within 5 per cent.

Page 31 The albedos calculated for this comparison were generated from the tracking of 1000 gamma-rays incident upon a 20 cm thick concrete slab for each incident energy angle combination.

The mean free path of a 2 MeV gamma-ray in concrete is about 10 cm. The fraction of reflected photons at a depth of 2 mean free paths is almost 100 percent (2). The concrete is taken to be that of ordinary concrete with a density of 2.35 g/cm3. The cutoff energy for history termination is 10 keV. The total albedos were found by numerically integrating the differential albedos over 5 exit polar and 5 exit azimuthal angle intervals.

Table 3: Total number albedo comparison 9 Study 2.0 MeV 1. 0 MeV 0.5 MeV 00 p 0.159 0.215 0.271 BR 0.162 0.221 0.268 R 0.164 0.207 0.275 60° p 0.317 0.351 0.420 BR 0.313 0.390 0.414 R 0.316 0.365 0.419 90° p 0.717 0.733 0.752 BR 0.724 0.744 0.734 R

Table 4: Total dose albedo comparison e Study 2.0 MeV 1. 0 MeV oo p 0.017 0.037 BR 0.020 0.040 CH 0.032 0.049 60° p 0.057 0.099 BR 0.055 0.099 CH 0.063 0.094 90° p 0.302 0.347 BR 0.303 0.355 CH 0.322 0.388 P Present paper -1000 histories BR Berger and Raso C9) R Raso (8) 0.5 MeV 0.071 0.074 0.076 0.143 0.146 0.136 0.387 0.395 0.456 Page 32 0.2 MeV 0.258 0.285 0.285 0.407 0.409 0.419 0.703 0.703 -----0.2 MeV 0.139 0.138 0.134 0.211 0.220 0.214 0.463 0.470 0.580 CH Chilton and Huddleston ClO> -numerically integrated Numbers unavailable Page 33 Use of Gamma-Rav Albedos in Sampling and Estimation Before any sampling is done a probability density function Cpdf) must be formed.

is given by a(Eo ,µ0 ,E,µ,¢>

The pdf formed from (19) where is the total number albedo defined by equation 3. For simplicity the parameters E., and µ:.;will be left off and p(E 0 ,fJc*E,µ,¢>

will be referred to as pCE,µ,¢>.

The definition of pCE,µ,¢> dE dfJ d¢ is then the probability that a particle with incident energy Eu and angle e 0 will be reflected into energy dE about E and into solid angle dfJ d¢ about e and ¢. However pCE,µ,¢> is a distribution function in three variables so we define the following pdf s p<E,Jl> = f T1 d¢ -T1 = 2 forrpCE,/J,fP>

dlfi pCE> = lo l p!E,/Jl d/J (20) ( 21>

Page 34 ( 22) We therefore sample for exit energy E' by inverting the equation {E' !1 = j O p CE> dE (23) Once E' is known then exit angle 8' or fl' = cose' is found from inverting t = 2 1:*p!E' .µ1 pCE') f, fl' = O pCfl/E' > dfl (24) After fl' is found along with E' the exit azimuthal angle ¢' is found by inverting 1¢* 2 0 p C E , , fl' , ¢ > d.f.* r = 3 pCE' ,fl') = 1¢' 2 O pc¢1µ',E'>

dtp <25)

Page 35 where r, , r2 and r3 are all independent uniformly distributed random numbers as defined earlier. The weight assigned to the particle is then ACE 0 ,µ0> or the total number albedo for a photon having incident energy E 0 and angle defined by µo = cos eo. This is the method used in the subroutine ALBSAM to sample for E, µ and ¢. However for purposes of program efficiency and computer memory the following is done. Letting Ej, µ k and ¢I be the maximum E, µ and ¢ in the energy or angular group designated by the indices j, k and 1 the following function is formed Since C<Ej,µk,¢1>

is actually a computed table the above is accomplished by j k = z L: L: j'=l k'=l where aj 1 k 1/ = value of diff. j = 1 to J the number of k = 1 to K the number of 1 = 1 to L the number of 1 L: l'=l a '1 L 1/1 I I II J number albedo at Ej, /Jk' exit energy groups; exit polar angle groups; azimuthal angle groups; ( 27) and ¢;I; Page 36 6E = constant exit energy interval; 6f1 = constant exit cose interval; 6¢ = constant exit azimuthal angle interval.

With this notation we can represent the total number albedo ACE 0 ,f1 0> as Therefore Cj, k,I I CJ 1 K,L is a cumulative distribution function (cdf> in three variables.

Thus it can be shown that the right hand sides of equations 23, 24 and 25 can be represented by !o"'i p(El dE = = Cj 1 k,L -Cj-1,k, L Cj)K 1 L -Cj-1) K 1 L Cj 1 1t,I -Cj-l)k 1 1 -Cj,k-1 1 1 + Cj-l,k-1, I Cj,k,L -Cj-1,k,L -Cj,k-1,L + Cj-1 1 k-1 1 L (28) (29) (30) To sample exit energy E from this table of C's we use equation 27 to find a j such that rl Cj,K,L/ cJ,K,L Dj ( 31)

Page 37 Linear interpolation then gives r, -Dj _I E = LlE + E. (32) I D* -Dj-1 I where D 0 is defined to be zero always and E 0 is the lower energy bound <EMIN> for the lowest energy group. Similar interpolations are carried out for µ and ¢ as well using equations 29 and 30. The program DOSE uses direction cosines to define the direction of the particle-track vector. Once the reflected direction given by µ = cosG and ¢ is sampled using ALBSAM it is necessary to translate the incident direction cosines ex, /.? and o into the exit direction cosines ex' , ;:;-' and o: This is accomplished in subroutine TRANS. Shown on figure 4 is the case where a particle is incident upon a surface with the Z-axis directed out from the surface and normal to it. Given are incident direction cosines O<, L? and polar angle of reflection e and azimuthal angle ¢. The emergent direction cosines are then given by ex' = sine cos< r;; + (J > /.j' = sine sin< 9' + (J > o * = cos e (33) (34) (35)

Page 38 z E y Figure 4: Direction cosines and angles of reflection where 'f is given by the relationships ex = cos;; cos/° o = -sin;° to obtain Page 39 (36) where S = 1 for /.] 0 and S = -1 for /.J < 0. Using figure 4 to visualize the opposite case where the Z-axis is directed into the surf ace the emergent direction cosines are found by the similarly derived expressions

()( I = sine cos (,9-' -¢) (37) tJ* = sine sin<JJ-¢>

(38) o' = -cose (39) where the angle 9 is again given by equation 36. In the case where the X-axis or Y-axis is normal to the reflecting surf ace the same relationships are equivalent expressions.

used to obtain To estimate the scattered dose at a point detector resulting from a gamma-ray incident upon a reflecting surf ace the photon must be weighted by the differential albedo and in the case of a number albedo an exit energy must be sampled. This is accomplished by the use of function subroutine DNALB.

Page 40 Given the incident direction cosines and those exiting in the direction of the point detector the subroutine UNTRAN, which does the reverse of TRANS, returns µ and ¢. The cosine of the exiting polar angle,µ, is found by setting it equal to the absolute value of the emergent direction cosine whose axis is normal to the reflecting surface. The angle;? is obtained by use of equation 36 and the azimuthal angle ¢ is found from either equation 33 or 37 depending upon whether the Z-axis is directed inward or outward from the reflecting surface as discussed earlier. Knowing E 0 , µ0 , µ and ¢ the indices m, i, j and k respectively are computed which symbolize the corresponding energy and angular groups. The indices m and represent E 0 and J.1 0 which as stated earlier are left off for purposes of simplicity.

The differential number albedo in solid angle is then found from = CJ,k 1 I -CJ)k-1)1 -CJ)k)l-1 + CJ)k-1) 1-1 2 ¥ fl¢ (40) Given the exiting direction defined by µ' and ¢' exit energy E' is sampled by inverting the equation

{E' Jo p(E,fl' ,¢'> dE' r = P<fl', ¢') {E' = }o pCE/fl' ,¢' > dE where the right hand side is represented by = Cj,k,I -Cj,k-1,I -Cj,k,l-1 + CjJk-IJl-1 cJ,1<,1 -cJ,k-1,1 -cJ,1<,1-1

+ cJ,k-1,1-1 Page 41 (41) (42) Exit energy E' is then interpolated the same way as described by equations 31 and 32 except that D;, is now equal to equation 42.

Page 42 A MONTE CARLO MODEL OF DELAYED FISSION-PRODUCT GAMMA-RAYS Introduction and General Discussion The computer model DOSE, written in FORTRAN IV, simulates delayed fission-product gamma-rays escaping from the NCSU PULSTAR reactor. Due to the complexities of the problem (e.g. three dimensional geometry and multiple scattering>

Monte Carlo is the method used. Although general purpose Monte Carlo photon transport codes are already available it was decided that a program specifically written for this application would execute much faster, take less memory, and provide a better understanding of the process. Because of the first two advantages the nuclear engineering departmental VAX-11/730 mini-computer was used to perform the computations.

It typically took about 5 hours5.787037e-5 days <br />0.00139 hours <br />8.267196e-6 weeks <br />1.9025e-6 months <br /> to obtain a complete dose map of the reactor bay using approximately 200 point detectors.

With 30,000 histories the standard deviations of almost all detectors are under 10 per cent. The program DOSE is general enough to allow the user to obtain dose rate information at as many as 300 different locations at once. Standard deviations for every detector are also computed.

or out of the Detectors can be placed anywhere inside reactor building.

This includes the biological shield, control room, reactor bay, reactor Page 43 platform, loading lock, roof and ground. All that is needed are the detector coordinates with respect to the axes shown in figure 1. All parameters necessary to describe the geometry of the reactor building are specified in subroutine BAY. To change the configuration only BAY has to be altered. The dose rate is a function of the individual detector location.

This requires the estimating of dose rate due to gamma-rays scattering in many different locations which in turn requires time. In order to minimize CPU time there are seven estimating options. These are !BAY, !TRANS, ICSCAT, ISSCAT, IDIREC, !SKY and IALD2. These are set either to 0 or 1, 0 meaning yes and 1 options are defined as follows; meaning no. for IBAY=O The the contributions to the dose from gamma-rays scattering off the reacter bay walls <including control room and loading dock, roof and floor) and the exterior of the biological shield are computed, otherwise they are not, for ITRANS=O the dose transmitted through or backscattered from the bay roof is computed using statistical estimation C2> instead of albedos, for ICSCAT=O the dose from scattered radiation within the core is computed, for ISSCAT=O the dose from radiation scattered within the biological shield is computed, for IDIREC=O the direct contribution to the dose is computed, for ISKY=O the dose due to scattering in the Page 44 atmosphere

<skyshine) is computed in subroutine SKYSHI, and for IALD2=0 the dose due to backscattering off concrete is computed using the C-H formula instead of the number albedo table. The latter is accomplished through the use of function subroutine ALD2. These options allow the computation of dose rates in every location using this single program. They also allow determination of the major sources of radiation since each can be computed individually.

The options apply to all detector locations.

The most important part of the model is the source itself. Small changes in the material composition of the PULSTAR core have a relatively large effect on the resulting dose rate. Table 1 gives the volume fractions of the material making up the core. Gamma-ray attenuation and scattering cross sections are taken from Storm and Israel (3). The core is homogenized and buildup is accounted for by following several scattering events inside the core. Buildup in the core generally accounts for 50 per cent of the dose on the bay floor. For each scattering event the probability of escaping the core is computed and the gamma-ray goes through a series of reflections within the biological shield. This is actually an oversized duct problem for which albedos are idealy suited. Using the albedo table the reflections continue until the gamma-ray escapes the shield or its Page 45 weight or energy become so low that it may be terminated.

Once the escaped photon history is completed, the tracking of the gamma-ray resumes at the point before escape from the core. Using subroutine FSCAT the gamma-ray is forced to scatter again within the core until it too is eventually terminated.

On the average 2 to 3 core scatters are followed per history. The biological shield is thick enough around the base and at the top to allow negligible contribution to the dose from transmission through the shield. It is therefore treated as an opaque structure.

However for detector locations at the top of the biological shield or on the bay roof the dose is due almost entirely to direct, core scattered, and shield scattered radiation escaping from the top opening of the shield. Hence the options IDIREC, ICSCAT and ISSCAT are all set to zero in this case. Eventually a gamma-ray escapes through the top of the biological shield and enters the bay. The most likely point of incidence is with the concrete bay roof. For a thick roof ClO.O cm) and detectors located within the reactor building the escape of gamma-rays through the bay ceiling is not treated. The photon is weighted by the total albedo and forced to reflect off the various concrete structures are met. until criteria for particle history termination Subroutine ALBSAM uses the cumulative Page 46 distribution functions formed from the table of differential number albedos to sample for exit direction and energy. For lower bay detector locations there is negligible contribution from radiation penetrating the shield and therefore only IBAY is set to zero since also has negligible effect. For thin roofs <the Lawrence Livermore benchmark case has a roof of 5/16 inch steel plate) there is a significant reduction in the dose rate due to gamma-rays escaping through the ceiling. The albedo concept is then overly conservative in the case of reflected dose and cannot be used at all for transmitted dose. Statistical estimation is used to estimate from the ceiling instead <IBAY and ITRANS=O>

meaning that the dose is estimated from each scatter within the ceiling and from each reflection within the rest of the bay. Although not allowed to escape the building the gamma-ray can still escape back into the bay as it can contribute further to the dose there. The same procedure is followed for detector locations on the roof where the transmitted dose must be computed.

If the detectors are located in such a way that contributions from direct. core scattered.

and shield scattered radiation are possible then IDIREC. ICSCAT. and ISSCAT must be set to zero as well.

Page 47 For detectors located outside the reactor building at ground level the dose is due almost entirely to backscattering in air. The NCSU reactor bay walls are constructed of concrete 12 inches thick which allows negligible transmission.

The concrete ceiling however is catacombed with the thinnest sections being approximately 3 inches thick. A conservative approximation is made in this case by treating the whole roof as if it were 3 inches thick. To estimate the dose the photon is tracked in the same manner as that used to compute the transmitted dose to detectors located on the roof with the exception that each time the gamma-ray is in the direction of escape through the ceiling the probability of escape is computed, a pathlength is sampled, and the gamma-ray is forced to undergo a series of scatters within an infinite atmosphere until termination criteria are met. For this case ISKY is set to zero which means that only scatters occuring in the atmosphere contribute directly to the dose. Ground effects are neglected.

As expected this case uses the greatest amount of CPU time per history. Once the gamma-ray scattering in the air is terminated the tracking resumes at the point just before its escape from the ceiling. Much has been said about termination criteria.

There are three ways in which the tracking of a photon in a particular region is stopped. The first way is to Page 48 terminate tracking of the photon if its energy falls below the cutoff energy of 20 Kev. When this cutoff energy was lowered to 10 KeV a negligible increase in the dose was observed.

The second way involves Russian roulette <2,19>. Each time the photon weight WC decreases it is compared to its preyious weight or the average weight WAVG before the decrease.

If WC is a certain fraction WMIN below WAVG then Russian roulette is used. A random number r is computed and if it is less than WC I CWAVG x WMIN> the photon weight is increased to WC = WC x WAVG x WMIN, otherwise it is terminated.

Function subroutine

!KEEP accomplishes this task. The optimal WMIN was found by trial and error to be approximately 0.05 for DOSE. The third method of termination is used during dose estimation.

If the maximum dose contribution to all detectors falls below 0.1 per cent of that already accumulated during the history, termination occurs and tracking resumes where it left off.

Page 49 Source Intensity and Energy Distribution The gamma activity is directly proportional to the dose rate. It is determined here by the expression given by Glastone <28) where A;. 0.7 p <t-0.2 _ t-0.2> 0 A = gamma activity in curies; P = average power level in watts; t = number of days after shutdown; tc = number of days after startup. (43) This equation is accurate up to a factor of 2 or less from 10 seconds to several weeks. Due to the fact that the PULSTAR has no fixed operating cycle the average power assumed is 1 MW even though in a typical month such as December 1983 the PULS TAR operated at the peak steady state power of 1 MW for only 44 hours5.092593e-4 days <br />0.0122 hours <br />7.275132e-5 weeks <br />1.6742e-5 months <br />. The time after startup is taken to be infinite.

Thus the reactor is assumed to have reached equillibrium power before shutdown.

These are the most conservative assumptions made in the analysis.

The actual fission-product activity could be much lower. The time after shutdown is taken to be 10 minutes. Substituting these values into equation 43 an activity of 1.9 x 10 6 ci is obtained.

from the loss of water Livermore 2 MW reactor was Page 50 As a comparison the activity experiment conducted at the estimated to be 3.0 x 106 ci after integration of the above equation over its power history. This is 20 per cent below the activity calculated for this reactor using the above assumptions.

The expression used for the approximation of the gamma-ray spectrum is taken from Maienschein (5). It takes the form N<E> = 7.4 exp(-1.10 E> I MeV C44> where N<E> is the yield and E is the gamma-ray energy. The constant corresponds to a total fission-product gamma-ray energy release of 5.9 +-0.7 MeV/fission for times after fission from 1 sec to 108 sec and gamma-ray energies above 2.8 MeV, plus 0.6 +-0.6 MeV/fission for energies below 0.28 MeV, plus 0.3 +-0.2 MeV/fission for times less than 1 sec. The exponent is the same as that used for the prompt fission gamma-rays.

Although roughly approximate we are only concerned with the exponent since the equation is to be formed into a pdf. Also this equation was compared to other equilibrium fission-product gamma-ray spectra <13) which provided energy release rates for seven energy groups along with the effective energy of the group. Converting this energy release data into the yield defined above Page 51 showed that equation 44 predicted a much harder gamma spectrum.

It is therefore thought that its use is at least a conservative approximation to the true fission-product gamma spectrum.

To sample from equation 44 we convert it to the pdf pCE) = 1.1 expC-1.1 E> ( 45) and invert its integral to arrive at the simple expression E = EMIN -ln 51 1.1 (46) The albedo table is good for energies up to 6.0 MeV. Energies sampled above this are set equal to 6.0 MeV. Photons having an energy above 6.0 MeV account for less than 0.2 per cent of the total sampled. The analytical model described later showed that the dose on the bay floor for purely ceiling scattered photons between the energies of 6.0 and 10.0 MeV contributed less than 1 per cent of the dose produced by those emitted between 0.01 and 10.0 MeV. In order to estimate the dose rate at times other than the reference time of 10 minutes after shutdown a log-log plot of data measured by Maienschein is provided in figure 5. The figure gives the fission-product gamma-ray activity versus time after shutdown for a reactor which operated for 1000 hours0.0116 days <br />0.278 hours <br />0.00165 weeks <br />3.805e-4 months <br /> at 1 watt. We are making the conservative Page 52 assumption that the fission product gamma-ray energy spectrum remains constant with increasing time after shutdown.

From a comparison of Maienschein's gamma-ray energy release rate to gamma-ray activity data, both as a function of time after shutdown, it is seen that the energy spectrum becomes softer with time. Table 5 below uses Maienschein's data to give the relative magnitude of the dose rate for different times after shutdown.

Table 5: Relative dose rate magnitudes for different times after shutdown Time after shutdown 100 sec 10 min 1 hr 1 day 1 wk 1 mo 1 yr Relative Dose Rate 1. 3 1. 0 0.73 0.29 0.14 0.059 0.003 Page 53 -I t-t-a: ;. -1010 -I u w lll ->-t--> -t-u a: 10 9 >-a: 0::: I a: a: L!l 10 8 10 2 20 3 20' io 5 20 6 io' TIME AFTER SHUTDOWN CSECl Figure 5: Total fission-product gamma-ray activity versus time after shutdown for reactor operating at 1 watt for 1000 hours0.0116 days <br />0.278 hours <br />0.00165 weeks <br />3.805e-4 months <br /> Importance Sampling From Source Angular Distributions The fission-product decay photons Page 54 Spatial are emitted isotropically.

However from figure 6 shown below it is seen th.p.t photons emitted in the solid angle subtended by the biological shield opening have a much better chance of escape than those incident upon the shield wall. It is then advantageous to put more importance on the sampling of photons in that direction.

Discrete importance sampling is used to accomplish this. Given the pdf p(x) defined over the interval [a,bJ we form a modified distribution p'(x) = c pCx> I(x) where c is constant and I(x) is the importance function.

The constant c is found from the property of pdf's defined by the equation l = lbp' <x> dx = c lbp<x> I<x> dx (47) Breaking the distribution p'Cx> into discrete intervals, each having importance In the above can be written as

-I STD. c 0,\, ,,::_ -.. ' <:;' .... . ' <l ' ...: !::.

  • BJ.,,!YYTES CONCRETE CONTROL ROD ACTUATORS TANK Z=O -6 , ' 4 ..> ,'; ' --ct_ Figure 6: Geometry of core and biological shield Page 55 Page 56 1 = c l I/ 1 p<xl dx + ... + In 1 p<xl dx x, 6Xn + .*. + IN l p<xl 6xN dx I N = c L: In l p<x> dx n=l 6Xn N l la Xn pCx> dx f.Xn-1 I = c L: In -a p<x> dx n=l N l Pn -Pn-11 = c L: In (48) n=l where 6xn = interval defining the nth region; In = importance value assigned to the nth region; N = total number of regions. The constant c is then given by 1 c = (49) N L In < Pn -P n-I> n=l To sample from p'<x> an mis found such that Then m-1 f, x r = c L I n C P n -P n -J ) + c Im p C x ) dx n=l x,,,_, (50)

Page 57 =ct:: In <Pn -Pn-/l + c Im <P<x> -Pm-/l (5ll Solving for PCx> we obtain P<x> 1 I m-1 ! = ---5-c} In <Pn -Pn-1> +Pm-I cim fi='1 (52) To obtain x one must invert the function P<x> given by P<xl =

dx ( 53) In this case x =µ, the cosine of the angle with respect to the Z-axis, a = -1, b = 1, and p<µ> = 112. The value µ is then given by p = 2 P

-1 (54) To account for the biasing the photon is weighted by w p<µ> = p' <µ> 1 = (55) c Im The rotational angle is interval [-fT,iTJ since their importance sampling here. Page 58 sampled uniformly on is little to gain the from This method may be used to importance sample source direction assuming the appropriate values of In are known. The importance of the individual directional regions can not be precisely determined until the problem is solved. The following procedure however will yield adequate estimates of these values. The dose and its absolute standard deviation is computed for one detector at each interval of A test case of a few histories C5000) is run and the standard deviations of the dose originating from source photons emitted in each interval is observed. Since the variance of the total dose is the sum of the variances of the individual doses it can be determined which interval is causing the greatest statistical error. This interval is also most likely to be one contributing significantly to the dose, hence it is an important region. The importance value In is then increased for that region. This method of trial and error may not be the most eloquent but with a little experimentation it works quite well. The standard deviation can be reduced by as much as 20 per cent for the same number of histories. The importance values In are read in as part of the input data into the array RIMPFCN>. Page 59 The decay fission-product gamma-ray source distribution is assumed to be identical to the power distribution of the reactor during steady state operation. Since the PULS TAR is a heterogeneous reactor the distribution is not entirely smooth. The horizontal distribution generally follows that of a cosine with peaks and dips due to the reflector, water channels and control rods between fuel elements. Fortunately it was found that the dose rate varied little as the horizontal power distribution was changed. It is therefore adequately approximated by the function p( x) = k cos ( ff x ) 2 <H + h> (56) where p<x> =horizontal power distribution in x <and y); k = constant; H = halfwidth of the core; h = extrapolated distance. The variable x is defined on the interval [-H,HJ and the extrapolated distance is calculated to be 6.7 cm. The dose rate is much more sensitive to the axial power distribution since photons originating in the top part of the core have a much better chance of escaping the shield. The equation used to describe this distribution is taken Page 60 from Tong and Weisman C27) and takes the form ,(rr k" z) (k" z) p<z> = k --cos 2 L L (57) where pCz> = axial power distribution; k', k" = constants; L = half height of core. This equation is a better approximation to an axial power distribution skewed by rod insertion. The variable z is defined on the interval [-L,LJ. The constant k" is found from the transcendental equation FCN/Z) = Cl.1585) k" I sink" axial peak power = (58) axial average power where FCN/Z) is the axial nuclear peaking factor and is equal to 1.51 for the PULSTAR (20). This gives a value of 1.227 for k". Shown on figure 7 is a plot of equation 57 superimposed on measured axial thermal flux data taken by a student during a laboratory exercise by activation of a copper wire. The equation predicts the thermal flux in the center quite closely but does not predict the rise in the flux at the ends of the core due to the reflector. This deviation however represents a relatively small area under the curve and is therefore neglected. -Xnld 18IX8 lD CT') o CT') I lD I O' I N -0 I U1 w :r: u z -z 0 -t--U1 0 CL ...J CI: -x CI: Figure 7: Relative axial thermal flux distribution Page 61 Page 62 Source positions sampled from the above distributions without biasing yield dose rates with very high variances and an unreasonable amount of machine time must be spent in reducing them to an acceptable level. This is due to the fact that low energy photons <which make up the majority of those sampled) have an average pathlength of only a few centimeters inside the core. Therefore many photons sampled towards the small probability of towards the surface. center of the core will have a very escape compared to those sampled To reduce the variance we wish to make each history contribute as equal a fraction to the dose as possible. Using z to represent x, y and z and p(z) to represent p(x), pCy> and pCz) we modify the pdf p(z) representing equations 56 or 57 to obtain the modified pdf p'(z) = p(z) exp<-pc<E> ID -zl> IN where /lc<E> = attenuation coefficient for the core; D = the z Cx or y> coordinate of the core surface the _photon is exiting; N =normalization constant causing .f p'(z) dz= 1. (59) Page 63 Equation 59 is basically the original pdf modified by the probability to exit the core unscattered. This modification is idealy suited to the unscattered gamma-rays emitted from the core but also works very well when the scattered component (buildup> is included. The weight to the gamma-ray is then W =WC x p(z) I p'(z) where WC is the photon's current weight. It was found that equation 59 had to be used with caution as its misuse could cause higher variances than that produced by the unmodified pdf. Thus the parameter D is found as follows. In order to determine which side of the core the source photon is exiting, its direction and position must be known. Therefore the source position is sampled from the unmodified pdf first. Once D and the direction are determined the same random number is used to sample from the modified pdf to obtain a coordinate which locates the photon emission closer to the surface it is exiting. Because p'Cz) is rather inefficient to sample from using a rejection technique another method is used. Subtracting the random number r from the cdf created from p'(z) one obtains g(z) = fz p'(z') dz' -r = 0 -L (60) Page 64 Brent's algorithm is then used to find the zero of g(z). This involves a combination of Newton's method along with bisection. There is no concern with multiple roots since g(z) has only one zero within the interval [-L,LJ. The same method is used to sample from equation 57. The function subroutine RFISM <Root Finding Importance Sampling Method) performes this task. Normally it takes from 3 to 4 iterations to reduce the relative error to less than 0.01 percent. Page 65 Geometrical Considerations Th.ere are many structures inside the reactor building that are not modeled in DOSE. Some structures, such as the control rod drive mechanisms located above the core, and the ceiling structure, could have a significant effect on the dose rate. Where approximations are made care is taken to keep the model on the conservative side. The control rod machinery for instance, will tend to attenuate many of the photons streaming out through the top of the biological shield and the structure of the bay ceiling will tend to trap photons directed upward and cause less reflection towards the control room and loading dock. Furthermore windows and wood doors are not modeled and scattering in air while in the reactor building is neglected. One aspect of the geometry which makes it much simpler to model is the fact that there are no curved surf aces and that almost every planar surf ace is perpendicular to one of the coordinate axes. The two exceptions to the latter are the base of the biological shield and the back of the control room. These are approximated by surfaces perpendicular to the axes. This makes it somewhat easier to compute angles of incidence and reflection. The direction that the wall is facing is kept track of by the variable !WALL corresponding to necessary part UNTRAN. which takes the 6 sides of the input Page 66 on the values of 1 through 6 of a cube. This is a to subroutines TRANS and All data names of structure locations start with a P for position and end in N, S, E, W, T or B for north, south, east, west, top and bottom. For example PBN stands for position bay north which indicates the position of the north bay wall with respect to the Y-axis. This data is all initialized in a block data subroutine. Modeling the reflecting surfaces is made simpler by the use of the two subroutines BOX and CUBE. Similar subroutines are used in MORSE-CG. Given a position and direction within a box or room, BOX returns the pathlength to intersection with the wall, the angle of incidence, and the direction with which the wall is facing CIWALL>. If the initial position is located outside of the box then the returned values pertain to the far side of the box. CUBE on the other hand computes the pathlength to intersection with a solid cube in space, the angle of incidence, the indicator !WALL, and the indicator IHIT which is 0 if the cube is hit and 1 otherwise. For both subroutines the caller supplies the 6 wall positions of the box or cube. Page 67 The biological shield itself is modeled quite accurately by a series of cubes surrounding a box representing the pool. The model is perhaps a little more detailed than necessary but it is done so that the effects of shadowing on the reflected gamma-rays can be observed. There are several modifications that could be made to DOSE and its subroutines to make it much more general (e.g. inclusion of curved surfaces>, but that is not its purpose. It is sufficiently general to allow easy modifications of the model and yet specific enough such that computing time is not wasted in decision making over various and unnecessary options. Page 68 Estimating the Dose As discussed in the Introduction General Discussion section there are a number of estimating options that can be used to include the contributions to the dose from different sources or compute them by another method. These options and how they are set for the different cases are explained there as well. Also discussed in the section called, Use of Gamma-Ray Albedos in Sampling and Estimation, are the equations used to form the number albedo distribution functions and their representation by discrete tabulated data. This section details the manner in which the dose rates and their associated statistical variances are obtained. There are three different methods used here to estimate the dose rate from scattered gamma-rays. The first method uses the number albedo table. This estimates the dose from a gamma-ray reflecting off a concrete structure and perhaps penetrating some barrier material. Once the appropriate differential number albedo DNALB and the sampled exit energy E are selected as described by equations 40, 41 and 42 the dose rate due to the reflection is computed by D = s WC DNALB exp(-µ( E> pcµ ,(n) K< E> I d2 ( 61> where S =source intensity (photons/sec>; WC= photon weight before reflection Ccm-2); DNALB = differential number albedo in solid angle; µ<E> = attenuation coefficient of barrier material Ccm-1>; pcµ,¢> = pathlength through barrier (cm); KCE> = flux to dose conversion factor; Page 69 lld 2 = solid angle subtended by point detector of unit area distance d away. Note that S WC is actually the current due to photons having incident angle e 0 on an elemental area of the reflecting surface. It has the units of photons per crn2 per sec. The approximation of l/d2 for the solid angle subtended by a unit area point detector is valid for large d but can cause problems if the point of reflection or scatter is very close to the detector. To prevent erroneous results all detectors are placed a few centimeters away from the scattering media. The second method of estimation uses the dose albedo as given by the C-H formula, equation 9. Since an exit energy is not known it is not possible to compute the attenuation coefficient. Therefore the transmitted dose through a barrier can not be computed. The dose rate is given by Page 70 (62) Note that S WC KCE 0) is the current incident on the reflecting surface in dose units. The third way of obtaining the dose rate is by a method known as statistical estimation. This technique is used to estimate from the scattered radiation within a scattering medium. It is given by (63) where pCE 0 ,9 5) is the probability per unit solid angle that a gamma-ray having incident energy E 0 will compton scatter through an angle 8 5* This is given by K<E::i, 8 5 > I oS or the differential over the total Klein-Nishina scattering cross section per electron. Supposing there are J reflections or scatters in this history then the dose rate computed for this history is X = t Dj < 64) J=l Now supposing there are N particle histories generated the total dose rate is given by u = 1 N Page 71 (65) where U is an estimate of the expected value of the dose rate. An estimate of the variance about the expected value is then given by (66) For each detector a sum of the contributions as well as a sum of the squares of the contributions from each history is stored to compute the variance. Page 72 BENCHMARKING A Comparison to Simple Analytical Solution As a preliminary benchmark for DOSE a simplified problem is formed, namely to compute the dose rate on the bay floor due to multiple scatters in the bay ceiling only. This excludes all secondary radiation exiting from the biological shield. Since almost all of the gamma-rays escaping from the shield unscattered will exit through the top of the core the horizontal source distribution will have little effect on the dose and is therefore neglected. Using equation 44 for the energy distribution and equation 57 for the axial source distribution we obtain the differential source term in E and z or photons S<E,z> = S p<E> p(z) (67) sec MeV cm where S is the source intensity in photons per second. If the rectangular shield opening is approximated by a circular opening of equal area then the gamma-rays shining on the ceiling is in the shape of the cone illustrated in figure 8. The differential dose incident upon elemental ceiling area dA is then d(z,E 0 ,8 0) where 1 = distance between source emission point to d.A; K<E> = flux-to-dose conversion factor; µc<E> = attenuation coefficient of the core; L = half height of core. Page 73 (68) The height, h, of the ceiling with respect to the center of the core is very large compared to the height, 2L, of the core; therefore, 1 can be approximated by and the variation in the angle of incidence eu can be neglected as the source location varies along the Z-axis. One is then able to analytically integrate the above equation over the height of the core to obtain (69) Multiplying dCE 0 ,8 0) by cose 0 we obtain the differential current in dose units incident upon dA. The differential dose at a point detector due to the reflection off of dA is given by I I I I h I g \ I \ m I \ ; ! I \ I I --4----...J d <POINT DETECTOR CORE \_REFLECTING AREA OF CEILING AXIAL SORCE DISTRIBUTION Figure 8: Geometry of simplified problem Page 74 Page 75 (70) where a 0<E 0 ,e 0 ,e,¢> is given by the C-H formula and m is the distance from dA to the detector. All that remains to be done is to integrate d<Ec,6 0 ,e,¢> over incident energy E 0 and the reflecting area N R2. The differential area dA as seen from figure 8 is given by dA = r d¢ dr. The total dose rate at the detector is then D = (EMAX {d+R 2 }o }d-R x d¢ dr dE <71> where the maximum azimuthal angle is a function of r and is given by -I cos (72) The cosine of the scattering angle is needed for tha C-H formula and is given by equation 5 stated here for convenience cos G 5 = sin 9 0 sin 9 cos i -cos .9 0 cos 9 ( 5) where in this problem cos e = g/m; sinG = rim; sin8 0 = s/l; cosec = h/l; and m = V r 2 + g2; s = Vd 2 + r2 -2 d r cos¢; 1 =vs2 + h 2. Page 76 The integration in equation 71 is carried out numerically using Simpson's rule. A Monte Carlo simulation of the same problem was made using both the C-H formula and the number albedo table. Below is a comparison table of some typical results as compared to the analytical solutions. Ideally the Monte Carlo run using the C-H formula should agree with the analytical solution exactly. They are indeed within the standard deviation computed. The results using the number albedo table are from 10 to 25 percent less than those of the other two methods. From table 3 it is seen that values for the dose albedo computed using the C-H formula are also generally higher than the literature. The inaccuracies of the C-H formula and proposed modifications to it have been studied in the past (24,26). Deviations in this Page 77 differential albedo with other data can be as great as 50 per cent or more. It is therefore not surprizing to find such variation in the results shown below. Page 78 Table 6: <mRem/hr) Dose rate comparison to analytical results Bay Control Reactor Method Floor Room Platform Monte Carlo Number .f\lbedo <1000 hist.) 22.3 +-3.6% 31. 8 +-4.1% 62.1 +-3.6% Monte Carlo Dose Albedo <C-H Formula, 1000 hist. ) 25.5 +-3.1% 41. 7 +-3.4% 77.5 +-3.3% Analytical Solution 25.9 42.7 78.9 Page 79 A Comparison to the Livermore Experiment As a final benchmark for DOSE a series of Monte Carlo runs were made those measured conducted at "Loss-of-Water to compare the calculated dose rates with during the loss of water experiment Livermore. Experiment This at the report, titled the Livermore Pool-Type Reactor," (4) gave dose rates measured at 10 minutes and 1 hour1.157407e-5 days <br />2.777778e-4 hours <br />1.653439e-6 weeks <br />3.805e-7 months <br /> after shutdown. Measurements were taken with both 0.0 and 1 foot of water over the core. The dose rates quoted here are the ones taken with 0.0 foot of water over the core. Although this experiment is similar to the type of accident modeled for the PULSTAR there are some major differences. First of all the Livermore reactor is rated at 2 MW and has a uranium fuel enrichment of 90 per cent. Unlike the PULSTAR this core is constructed of plate type fuel elements of uranium aluminimum alloy. The core and fuel element dimensions obtained from the Livermore reactor design data <12) yield weight percents of the materials comprising the core to be 79 percent aluminum, 18 percent water, and 3 percent uranium. There are also significant differences in the geometry of the biological shield and building. As shown in figure 9 the reactor building .,._ _____ sc' dia. BA 'I ROOF ACTUATORS J 11 5 TEEL (j) 3 PLATFORM __ .... 91_6.._'_.7 .... ,,---__, 3 EI (... dia. -I BIO. 5 HIE L D 26 1 6 II -( *-f_-BA Y FL C OR 4 core -0 . --t> *. 60' 90' Page 80 \:3 I 2C ' (§} -------15C' --------"'L-" Figure 9: Lawrence Livermore reactor building Page 81 is a domed cylindrical structure constructed of steel plate. The biological shield is also cylindrical in shape. In order to model the Livermore reactor facility using DOSE the reactor building is approximated by a rectangular steel structure of the same height and floor area. The cylindrical reactor pool is approximated by a rectangular pool of equal height and volume. Because thermal flux distributions for this core were unavailable the shape is assumed to be the same as that used for the PULSTAR. The fission product activity is taken to be 3.0 x 106 ci which was estimated by those conducting the experiment. Five runs are made using the above activity as the source intensity. Both the calculated and measured values are listed in table 7. To identify the representative locations of the computed and measured dose rates table 7 lists numbers under the label "location" which correspond to the circled numbers shown on figure 9. The best agreement is obtained on the reactor balcony which is not in direct line-of-sight with the core. Here the values are within 13 per cent of each other. The dose rates on the roof compare somewhat less favorably. The maximum computed dose rate of 290 +-12 mRem/hr is more than twice the measured value of 105 mRem/hr. The Page 82 Table 7: Dose rate comparison to measured values Measured Location Dose Rate 1 (Reactor Platform) 195 mRem/hr 2a (Bay Roof) 105 mRem/hr 2b 46.5 mRem/hr 2c 0.58 mRem/hr 2d 0.84 mRem/hr 3 (Line-of-Sight) i 4750 Rem/hr ii 2500 rad/hr 4 (Bay Floor) iii 370 mRem/hr 5a (Outside) 80-60 mRem/hr 5b 55-40 mRem/hr 5c 35-30 mRem/hr 5d < 30 mRem/hr i Estimated value iv v Computed Dose Rate 220 +-6% mRem/hr 290 +-4% mRem/hr 24 +-2% mRem/hr 0.27 +-7% mRem/hr 0.67 +-8% mRem/hr 2925 +-4% Rem/hr 2135 +-4% rad/hr 70 +-6% mRem/hr 720 +-8% mRem/hr 72.3 +-18% mRem/hr 33.2 +-17% mRem/hr 18.7 +-15% mRem/hr 12.8 +-15% mRem/hr ii Measured 1 hour1.157407e-5 days <br />2.777778e-4 hours <br />1.653439e-6 weeks <br />3.805e-7 months <br /> after loss of water iii Maximum dose rate measured on bay floor iv Unconservative case (allows for escape through roof> v Conservative case (no gamma-rays escape bay> Page 83 taking place is probably due to some added attenuation in unmodeled structures such as control rod actuators, structural beams, etc. The calculated dose rate line-of-sight of the core is 20 feet above and in 38 per cent below the 4750 rem/hr estimated in the experiment. This apparent underestimation needs to be looked at in light of the fact that the quoted value is estimated and not actually measured. The dose rate measured at this same point 1 hour1.157407e-5 days <br />2.777778e-4 hours <br />1.653439e-6 weeks <br />3.805e-7 months <br /> after shutdown was 2500 rad/hr. If we use table 5 to account for the decrease in activity one obtains for the computed value 2135 +-85 rad/hr which is now only 15 per cent below the measured dose rate. This indicates that the difference might result more from an error in their estimation rather than an error in the calculated value. The dose rates on the bay floor are computed in two ways. The first allows for the fact that many of the gamma-rays escape through the relatively thin steel roof. The second and more conservative way is to assume that there are more structures in the* bay preventing escape than those modeled. The control rod actuators for instance are located directly above the core and will tend to backscatter many of the gamma-rays escaping the biological shield. Thus a conservative value is obtained by that none of the photons escape the building. assuming The true Page 84 value should lie somewhere in between as evidenced in table 7. The PULSTAR on the other hand has very thick bay walls Cl ft) and the ceiling is considered to be thick also in the case of detector locations within the reactor building. The computed dose rates should then be conservative in the case of the PULSTAR. As a final comparison dose rates are computed 3 feet above ground level outside of the reactor building. The effects taken into account are transmission through the bay walls and skyshine. The detectors are placed from 60 to 150 feet out from the center of the reactor building. The computed dose rates decrease from 70 to 10 mRem/hr with standard deviations of approximately 15 per cent. The dose rates measured over the same distances decreased from 80 to 30 mRem/hr. The calculated dose rates would probably be much higher if the shape of the actual domed roof were taken into account as some contribution from transmission through the ceiling would result. It must be kept in mind when comparing these values that several assumptions and approximations have been made in this simulation of the Livermore experiment, some of which are not made in the PULSTAR model. The calculated dose rates however generally agree with the measured values to better than 50 per cent. It is therefore concluded that Page 85 the model is indeed a good one and should give even more accurate results for the PULSTAR for which it was written. Page 86 RESULTS AND CONCLUSIONS The introduction to this paper pointed out the need to answer some basic questions concerning the dose. They basically ask how much, how long, and where. The analysis is made. keeping the above questions in mind. Detectors are located such that isodose lines, lines of a constant dose rate, may be established. Approximately one point detector is used per square meter of analyzed surface. All units are in either mRem/hr or Rem/hr as indicated. A standard deviation is computed for each detector as described by equation 66. Those shown in the following figures however are obtained by averaging the individual deviations of the detectors within the area defined by the isodose lines. The first area of interest is the bay floor. The elevation of all detectors is z = 0 or 0.61 meters above the floor. Figure 10 shows the dose rates computed using 30,000 histories and 5 hours5.787037e-5 days <br />0.00139 hours <br />8.267196e-6 weeks <br />1.9025e-6 months <br /> of CPU time. As seen from the figure the hottest spots <150 to 175 +-7% mRem/hr> are located east and south of the reactor. This is due to the off-center location of the core within the biological shield and greater distance between the reactor and south bay wall than between the reactor and north wall. The coolest <least exposed> area is oddly enough adjacent to the biological shield. Due to the fact that 75 per cent of Page 87 the dose is due to shine off the ceiling with the other 25 percent originating from the rest of the bay the biological shield shadows the area around it. Thus the north bay door provides the means of access to the bay floor with the least exposure to gamma radiation. As* a comparison the same analysis is made using the C-H formula as an estimator instead of the number albedo. The results are similar although about 15 percent higher on the average. The hot spots, as shown on figure 11, are predicted to be in the same locations and the same conclusion as that made above is drawn. One might note that the computed errors are slightly lower in this case since there is not the added statistical variance caused by sampling exit energy. Another point of interest is rates predicted on the bay floor. made concerning dose The average background on the floor while the reactor is operating is measured to be approximately 0.2 mRem/hr. According to figure 10 and table 5 it will take approximately 1 year before the maximum computed dose rate of 175 +-12 mRem/hr is reduced to a value of twice the present background. The dose rates on the reactor platform shown in figure 12 are computed for detectors located 9 meters above the bay floor. The magnitudes as expected are much higher due Page 88 to the reduced distance between the platform and bay ceiling. The elevation of the detectors are approximately head height if one is standing on the platform. Because much of the interior surf ace of the biological shield is visible from the platform there is a significant contribution to the dose rate due to gamma-ray reflections within the shield. This accounts for the rapidly increasing dose rates as the detectors are placed closer to the edge of the pool. Other than a location in direct line-of-sight with the core the platform offers the area of greatest exposure to the gamma field and is on the order of several Rem/hr. Figure 13 illustrates the expected dose rates found in the control room. The gamma-rays stream in through the two control room windows and door much as the sun shines in a window. The effect is to cause dose rates on the order of 250 +-13 mRem/hr next to the windows and door decreasing rapidly to a minimum of 5 +-2 mRem/hr in the far corner. The detectors are placed 1 meter above the control room floor. Although very high the dose rates will undoubtedly be reduced drastically by the placement of temporary shields (e.g. lead bricks) over the windows and door. Page 89 A similar effect is seen in the loading dock area as illustrated in figure 14. The projection is a small landing 12 feet above the bay floor and on the same level as the loading dock. The equipment room is located underneath. The loading dock is a non-radiation area frequently used as a passageway between buildings located north and south. Dose rates as high as 175 +-9 mRem/hr tapering off to a fraction of a mRem/hr are encountered in this area. As indicated previously attenuation through the door opening into the bay is assumed negligible. Other areas of interest and ones that may have to be reentered are the offices of the Health Physicist and Nuclear Reactor Operations located adjacent to and north of the reactor building. Due to the attenuation provided by the solid 1 foot thick concrete bay wall the dose rates predicted are on the order of 10-6 +-10% Rem/hr. This is comparable to the present background and is considered negligible. Dose rates found on the bay roof range from a few Rem/hr when not in line-of-sight of the core to 25 +-2 Rem/hr when directly overhead 16 meters above its midplane. Values are shown on figure 15. The dose rate is also computed at a location 7 meters above the midplane of the core and in its direct line-of-sight. This puts the detector on the same level as the top of the biological shield. A value of 230 +-10 Rem/hr is calculated. Page 90 The final yet very important area considered is outside the reactor building. As seen in figure 16 the isodose lines form concentric circles around the reactor building with the reactor core at the center. These dose rates are computed using dry air as the scattering medium. Ground effects are neglected. The dose is strictly due to skyshine as the bay walls are too thick to allow significant transmission. As illustrated the dose rate assumes a high value of 4 +-0.3 mRem/hr immediately adjacent to the building to a low value of less than 0.3 +-0.02 mRem/hr for distances greater than 80 meters from the core. Assuming that the natural background outside the reactor building is 200 mRem/yr Cthe average natural background for the United States> the dose rate 10 minutes after the accident would be roughly 200 times the background level at a location next to the building. In conclusion it can be said that the maximum dose rate occurs on the reactor platform. Depending upon how close one gets the dose rate can range from a fraction to several hundred rem/hr when located directly over the core. The bay roof is also an area that is subject to a high intensity gamma field, but this is not an area normally occupied by personnel. From calculations based upon the loss of water experiment conducted at the Livermore Page 91 Pool-Type Reactor the PULSTAR Safety Analysis Report <20) estimates the radiation level outside the control room window to be approximately 100 mR/hr with the pool drained. The corresponding exposure calculated here is 350 +-17 mR/hr or 3 and 1/2 times as high. Both calculations assume an power of 1 MW before shutdown. It must be understood however, that the values reported here are of necessity somewhat conservative. As stated earlier the fission-product gamma spectra used in this study is somewhat harder than that measured by Maienschein (5). From a quick glance into the reactor bay one will also notice many other attenuating structures Ceg. stairways, actuators, structural beams) besides the biological shield and bay walls. These are left out in the interest of a conservative analysis. As far as dose rates outside of the reactor building are concerned it should be noted that personnel located inside the building away from bay openings may be less exposed than those standing immediately outside due to the shielding provided by the building against gamma-rays scattering in air above the reactor building. NORTH--100 -125 a 8 % 125 -150:J:. 6 % SCALE 0 I 2 3 4 5 meters NORTH BAY DOOR 50-75:t 10% 25-50:t 10% Page 92 Figure 10: Computed dose rates on bay floor CmRem/hr) Page 93 NORTH 100-125 125-150+/-7 % 150-17S:t. 6 96 75-100 :t. 7 % 50-75 :t. 7 96 25 -50.:t. 8 96 SC AL E ------o 2 3 4 5meters Figure 11: Computed dose rates on bay floor using the C-H formula CmRem/hr) NORTH ,_-----..._ .,,, ,.,,,,. .,,,, .................... ' / ' , ' / .,...,---- ..... ,, ' ,, ,. ... ' I / ' \. .,,. l-10z5% , ' .,,. \. I / I I I I I l \ \ \ ' ' \ ' ' ' .. ' ' .. / ... 0.5: l.t.-5 % ,, --0.3 -0.5 :L 5 % SCA L E I ,, / / I ' I I I I / / / 0 I 2 3 4 Smefers Page 94 10-50.z.496 50-230 .t. 4 % Figure 12: <Rem/hr) Computed dose rates on reactor platform ....... w .. () 0 51 "O c rt ro p. p. 0 t:ll ro "'! Pl rt ro t:ll ..... ::s n 0 ::s rt "'! 0 ....... "'! 0 0 51 DOOR ---NORTH SCALE 0 I 2 CONTROL ROOM WINDOWS 3 ml'fer s Page 96 t.O ... t.O t.O ... " c c C\l I t.O ...... ---C\l -:x: ..... co Q: -c C) t.O

  • l Lw -...J \.) U) Figure 14: Computed dose rates in loading dock <mRem/hr>

Page 97 NORTH ) < I :t 10 % SC ALE 0 2 3 4 5mefESrs Figure 15: Computed dose rates on bay roof <Rem/hr> Page 98 1-2+/- 7% 0.5-1 t 7% EAST BROUGHTON OR. 10 20 Figure 16: Computed dose rates outside reactor building <mRem/hr) Page 99 LIST OF REFERENCES

1. L. L. Carter and E. D. Cashell, Particle Transport Simulation With the Monte Carlo Method, Technical Information Center, Oak Ridge, Tennessee (1975). 2. N. M. Schaeffer, Reactor Shielding for Nuclear Engineers, U. S. Atomic Energy Commission (1973). 3. E. Storm and H. I. Israel, Photon Cross Sections from 0.001 to 100 MeV for Elements 1 Through 100, University of California, Los Alamos, New Mexico (1967). 4. M. Knezevich, R. L. Kathren, 0. K. Helferich and K. R. Kase, Loss-of-Water Experiment at the Livermore Pool-Type Reactor, Health Physics Journal, 11: pp. 481-487 (1965). 5. F. C. Maienschein, Fission-Product Gamma Rays, Engr. 6. Compendium on Radiation Shielding, 1:76-84, Weslay, New York Cl968}. C. J. Everett and E. the Inverse of Distribution, USAEC Scientific Laboratory D. Cashwell, Approximation for the Klein-Nishina Probability Report LA-4448, Los Alamos (1970). 7. C. J. Everett, E. D. Cashwell and G. D. Turner, On a New Method of Sampling the Klein-Nishina Probability Distributions for All Incident Photon Energies Above 1 KeV, USAEC Report LA-4663, Los Alamos Scientific Laboratory (1971). 8. Dominic J. Reflection Nuc. Sci. Raso, Monte Carlo Calculations on the and Transmission of Scattered Gamma Rays, Eng., 17: pp. 411-418 (1963). 9. Martin J. Berger and Dominic J. Raso, Monte Carlo Calculations of Gamma-Ray Backscattering, Radiat. Res., 12: pp. 20-37 (1960). 10. A. B. Chilton and C. M. Huddleston, A Semiempirical Formula for Differential Dose Albedo for Gamma Rays on Concrete, Nuc. Sci. Eng., 17: pp. 419-424 (1963). 11. R. N. Blomquist and E. M. Gelbard, An Assessment of Existing Klein-Nishina Monte Carlo Sampling Methods, Nuc. Sci. Eng., 83: pp. 380-384 (1983). 12. J. B. Radcliffe, Jr. and E. E. Hill, Lawrence Radiation Laboratory (Livermore), Report UCRL-4919 Rev. (1960).

Page 100 13. Reactor Physics Co:nsta_nts, Second Edition, pp. 629-630 (1963). -. . I ---14. D. George, T. R. England, R. J. Labauve and c. W. Maynard, Photon Sources for Shielding Applications, Trans. Arner. Nuc. Soc., 35: pp. 463-464 (1980}. 15. E. M. Fournie and A. B. *Chilton, _Gamma-Ray Buildup Factors for Concrete Slab Shields Under Slant-Incidence Conditions, Nuc. Sci. Eng:., 76: pp. 66-69 ( 1980). 16. Kazuo Shin, Ryuj i Murakami"; Hiroaki -Taniut:hi, _ 'rumonori Hyodo *:and Yoshiaki'

. Oka, :;Measurements*-

of -Neutron and Garnma:-Ray Streaming _ in.i *:.:a Ca vi ty:::Duct x: System. -and ran Analysis. by* : an Albedo Monte Carlo Method, Nuc. Sci. Eng., 81: pp. 161-171 (1982). 17. Ivan Lux, A Handy Method for Approximate-optimization of Splitting in Monte Carlo, Nuc. Sci. Eng., 83: pp. 198-205 (1983). 18. Ainericari National Standa:rd Guidelines

  • on the -Nuclear Analysis and Design of Concrete Radiation Shielding for Nuclear Power' *Plants, American Nuclear Society, ANSI/ANS-6.4 Cl977}. 19. American:

National Standard _for Calculation and _Measurement of Dire<:t and Scatte_red Gamma Radiation from LWR Nuclear Plants, American*

Nuclear '.-Society., ANSI/ANS-6.6_.l (1979}. NCSV PULSTAR Safety Analysis, C_opy -*Number 80, of *Nuclear-
  • NCSU, -Raleigh, North Carolina*

(1967}. 7:J... M* . and D._ .J !_ ** ._ 0£trnnia _ <1960). , . . -J ... ,,: .. ' -... .' '!' 22. C. M*. ,_ Davisson and L. D. Beach, A Monte Carlo Study bf Gamma.Radiation,-Trans. Amer.

  • 391-(1962->.
23. A. B. Chilton, C. M. Davisson, and L. A. Beach, Parameters

_-.for .c

  • the C-H .. Albedo Formula for Gamma Rays Reflected from Water-, Concrete, Iron, and_ Lead, Trans. Arner. Nuc. _Soc. , 8: pp. 656-657 C 1965)
  • 24. A. B. Chilton, A Modified Formula for. Differential Exposure Albedo for Gamma Rays.Reflected from Concrete, Nuc. Sci. Eng., 27: p. 481 (1967).

Page 101 25. L. G. Haggmark, T. H. Jones, N. E. Scofield and W. J. Gurney, Differential Dose-Rate Measurements of Backscattered Gamma-Rays from Concrete, Aluminum and Steel, Nuc. Sci. Eng.,23: pp. 138-149 (1965). 26. M. B. Wells, Differential Dose Albedos for Calculations of Gamma-Ray Reflection from Concrete, USAEC Report RRA-T56, Radiation Research Associates, Inc. (1964). 27. L. S. Tong and J. Weisman, Thermal Analysis of Pressurized Water Reactors, Amer. Nuc. Soc. Pub., La Grange Park, Illinois, p. 255 (1979). 28. S. Glastone, Principles of Nuclear Reactor Engineering, Van Nostrand, Princeton, p. 120 <1955). 29. Wei-Chien Tung, A Guide to General-Purpose Monte Carlo Gamma-Ray Transport Code With INER-0409/MN-80 (1981). Computer Code MORSE-CG: Multigroup Neutron and Combinatorial Geometry, 30. P. A. Robinson, Jr., The Development and Application of a Coupled Monte Carlo Neutron-Photon Transport Code, USAEC Report UCRL-51234, University of California Lawrence Livermore Laboratory (1972). 31. C. D. Zerby, J. Agresta, et al., PHOTRAN, A General Purpose Photon Transport Program in Complex Geometry, Technical Report AFWL-TR-65171 <Vols. I-IV>, Union Carbide Corporation, Research Institute (1966-1968).

32. E. D. Cashwell, J. R. Neergaard, C. J. Everett, R. G. Schrandt, W. M. Taylor and G. D. Turner, Monte Carlo Photon Codes: MCG and MCP, USAEC Report LA-5157-MS, Los Alamos Scientific Laboratory (1973). 33. E. D. Cashwell, J. R. Neergaard, W. H *. Taylor, G. 34. D. Turner, MCN: A Neutron Monte Carlo Code, USAEC Report LA-4751, Los Alamos Scientific Laboratory (1972). J.W. Kimlinger, F. Plechaty, and SORS Monte Carlo Photon-Transport Code USAEC Report UCRL-50358, University Lawrence Radiation Laboratory (1967). J. R. Terrall, for the CDC 6600, of California
35. E .. D. Cashwell, C. J. Everett, and W. M. Taylor, A General Monte Carlo Photon Code, in A Review of the Monte Carlo Method for Radiation Transport Calculations
Retrieved from "https://kanterella.com/w/index.php?title=ML18312A310&oldid=1210051"