A Monte Carlo Radiation Transport

Example Calculation Using PFC

Radiation transport is fundamental to the analysis of the radiation environments produced by nuclear materials and processes. It is used, inter alia, for calculating the effectiveness of biological shields, tissue dose profiles associated with radiation oncology treatments, radiation environments in the transport of radioactive materials in both normal and accident conditions, critical mass and power distributions in fissile assemblies, self-shielding of fertile materials undergoing activation for the production of medical and industrial isotopes, and scattering environments produced in experiments with particle accelerators. The Monte Carlo method is a powerful tool for solving radiation transport problems of all kinds, but is particularly important for the solution of problems involving particle motion in complex geometries. Familiarity with the Monte Carlo method and its application to the Boltzmann transport equation is important for practitioners in the field of radiation transport. A Monte Carlo Primer, A Practical Approach to Radiation Transport (the Primer), and its companion A Monte Carlo Primer Volume 2, Solutions to Exercises (the Solutions), provide a basic description of the Monte Carlo method as applied to radiation transport plus an in-depth explanation of techniques commonly used to exploit the method in realistic situations.

The Probabilistic Framework Code (PFC), a Monte Carlo transport code written in the Fortran language, is available as a free download elsewhere on this web site. For a detailed description of this code, and examples of its application to a variety of radiation transport problems, see Chapters 5-9 of both the Primer and the Solutions. The following example illustrates the power of the Monte Carlo method for solving radiation transport problems, and the use of PFC for obtaining Monte Carlo transport solutions quickly and easily.

Consider a point, isotropic source of neutral, monoenergetic particles (which we can consider for the moment to be neutrons) located at the center of a sphere of pure scattering material that is ten mean free paths thick to the source particles. Assume that the sphere is surrounded by vacuum. Assume further that inside the sphere the particles scatter isotropically in the laboratory coordinate system and that their mean free path is constant. We wish to determine the angular dependence of the current of particles escaping from the sphere. For a discussion of the mean free path of a particle, isotropic scattering in the laboratory system, the vacuum boundary condition, and particle current see Chapter 3 of the Primer or see other references on radiation transport.

Let us assume the source is located at the origin of a Cartesian coordinate system and that the mean free path of the source particles is 1 cm. Because the sphere is a pure scattering material we know that all of the source particles will escape. Furthermore, because the sphere is relatively thick to the source radiation, and because the particles scatter isotropically within the sphere, the particle flux at points far from both the source and the outer boundary should be nearly isotropic. We expect the vacuum boundary condition to perturb this isotropy, however, and thus we expect the escaping current of particles to be anisotropic. To understand the following problem description we recommend the reader refer to Chapter 5 of the Primer. As is often the case with Monte Carlo programs, it is necessary to make changes to the PFC subroutines in order to solve specific problems. The subroutines being modified are described in the Primer.

In order to use PFC to score the current of particles leaking from the sphere of scattering material we will modify subroutine ‘Bdrx’ to score the leaking particles as a function of the cosine of the angle between the particle trajectory, given by the unit vector , and the outer normal to the sphere at the point of escape, given by the unit vector n. This geometry is described on page 199 of the Primer. The dot product p· is given by eqn 3.48 of that volume, where p = rn and r is the radius of the sphere in the present calculation. Thus to normalize the result shown in eqn 3.48 for this calculation the dot product must be divided by the radius of the sphere. We define ten equal intervals of the cosine over the range [0,1] and score the leaking particles in these bins. We must add common ‘Stats’ and the associated double precision declarations (see subroutine ‘Stats’) to subroutine ‘Bdrx’ as given in Table 5.19 of the Primer. Then, following statement number 15 of ‘Bdrx’ we add the following (or equivalent) lines of Fortran:

i=INT((x*u+y*v+z*w)+1)

bscore(i) = bscore(i)+1.

bsumsq(i) = bsumsq(i)+1.

In addition to these changes in subroutine ‘Bdrx’ we need to make a few changes in subroutine ‘Stats.’ In statements 22 and 24 in the PFC library version of ‘Stats’ given in Table 5.9 of the Primer we replace the subscript ‘1’ with an ‘i.’ Prior to statement 22 we add

DO i = 1,10

We then delete statements 25 through 27 and replace them with

END DO

Using the geometry and cross section definitions shown in Tables 5.23 and 5.24 of the Primer we are now ready to solve the present problem. Running start particles with a start random number of one gives the results shown in the table below. It is clear from these results that the leakage is a maximum at a cosine of one and a minimum at a cosine of zero; i.e., that the maximum in the leakage current per unit cosine interval occurs in a direction normal to the surface of the sphere. This is clearly shown in the normalized polar plot of similar results obtained using a large number angular intervals that is presented in the figure below. In this plot the radius r of the curve, , represents the normalized probability of a particle leaking from the sphere of scattering material at the angle , where is the angle between the particle trajectory and the outer normal to the sphere at the point of escape. The radial values are plotted at the center of the associated cosine intervals. The leakage probability is normalized so that the maximum probability is one. Because the uncertainties are relatively small in these results such uncertainties are not indicated in the figure, although some statistical fluctuations are apparent. We see that, not only is the leakage probability a maximum for directions normal to the sphere, but the probability decreases rapidly as the angle to the normal increases.

Let us suppose one could view this sphere from a distance with an imaging sensor that indicates the magnitude of the flux of particles striking each pixel in the sensor. Because the sphere is in a vacuum, the sensor would measure the number of particles leaking through the surface of the sphere as a function of position on the two-dimensional projection of the sphere perpendicular to the line of sight between the sphere and the sensor. This two-dimensional projection is a disk. The fluence measured at or near the center of the disk consists of particles that escape from the sphere in a direction almost normal to the sphere, while that measured at the outer edge of the disk consists of particles that escape from the sphere in a direction nearly tangent to the sphere. Therefore, based on the present results the apparent brightness of the disk would be greatest at the center and would decrease with distance from the center of the disk to a minimum at the edge. Although the analogy is not exact, when the disk is the sun, the escaping particles are visible photons, and the sensor is either a photographic plate or the human eye, this phenomenon is known as limb darkening.

Upper Cosine Boundary	Leakage Probability	Standard Deviation
0.1	0.003021	5.5E-5
0.2	0.01385	0.00011
0.3	0.02812	0.00017
0.4	0.04553	0.00021
0.5	0.06751	0.00025
0.6	0.09350	0.00029
0.7	0.12442	0.00033
0.8	0.16074	0.00037
0.9	0.20436	0.00040
1.0	0.25891	0.00044

This problem is similar to some of those discussed in the Primer and the Solutions. As is done in those volumes, let us consider some of the ramifications of the present calculation. In the scoring done here, any particles escaping from the sphere without experiencing a collision are combined with the scattered particles in the most forward-directed angular interval. However, the uncollided particles are exactly radially directed in this calculation and thus could be expressed as a delta function at in the polar coordinates of the figure. The curve of the angular dependence should thus properly have a narrow peak normal to the sphere. The area under this peak will depend on the fraction of escaping particles that are uncollided. In the present calculation, for a sphere ten mean free paths in radius, this fraction will be small.

If we vary the radius of the sphere we expect the shape of the angular leakage to change. As the sphere becomes smaller, the fraction of escaping particles that are uncollided will increase. Thus the narrow peak in the angular dependence of the leaking particles at will increase relative to the scattered contribution. In addition, the average number of collisions experienced by an escaping particle will decrease and thus the anisotropy of the scattered contribution will become increasingly forward peaked. In the opposite limit, as the radius of the sphere becomes much greater than the mean free path of the particles, the problem can be approximated by an infinite half-space of scattering material. This approximation is commonly made for calculating the angular dependence and polarization of photons escaping from stars.

In a physically realistic configuration a source of particles can never be a point. Instead, the source will have a non-zero radius and the uncollided contribution to the particle leakage in a configuration such as that modeled here cannot be a delta function in angle. For the present calculation the angular dependence of the escaping particles will vary not only with the radius of the sphere but also with the radius of the source. In practice the source may vary in strength per unit volume with distance from the center of the sphere.

All of these variations can be easily incorporated into the Monte Carlo model used in this example calculation, and the characteristics of the angular leakage can be examined as a function of the problem parameters. However, because all Monte Carlo results contain statistical uncertainties (over and above any systematic uncertainties) comparisons between the results obtained from similar Monte Carlo calculations can be misleading. In many cases it is useful to make use of the technique of correlated sampling to isolate and quantify the differences between such results. Correlated sampling is discussed in Chapter 9 of the Primer.

Return to Home Page.