Transport theory and diffusion theory are based on equations that describe the average behaviour of large numbers of neutrons as they interact with large numbers of nuclei. As the preceding paragraphs show the equations are complex and their solution by numerical means even more so. This is in contrast to the actual behaviour of individual neutrons, which is quite simple: they travel in straight lines between nuclei and they interact with the nuclei in a limited number of ways the probabilities of which are known. The transport equation has the effect of turning a large number of simple problems into a single complex problem. The alternative, which is to solve the simple problems, is the “Monte Carlo” method (Brown, 2012).

The basis of the method is to track neutrons one by one. The outcome of each event in the history of each neutron is chosen at random in accordance with the known probabilities. This is done for large numbers of different neutron histories, so large a number that when they are all put together they make up an estimate of the actual neutron distribution.

The method depends on a random number generator that produces a sequence of numbers Rn distributed uniformly at random in the range (0,1). For a “source-type” calculation (for example a subcritical reactor driven by a neutron source) the procedure is straightforward. A neutron is assumed to emerge from the source travelling in a direction (в, f) chosen at random (which is ensured by setting f = 4nR1 and в = cos-1 (1 — 2R2)) and with an energy E given by R3 = jf s(E’)dE’, s being the spectrum of source neutrons normalised to 1. The reactor is assumed to be made up of discrete regions (fuel, coolant, structure, etc.) each of which has a uniform composition and uniform nuclear properties. The neutron is assumed to travel in its initial direction for a number m = — ln R4 of mean free paths to its first interaction with a nucleus. In its course it may well cross one or more boundaries from region to region. The actual distance travelled is given by m = x1E(1 + x2£(2 + … Хі^і where x1 is the distance travelled in the region in which it was born, x2 the distance through the next region, and so on, and xi is the distance travelled in the region i in which it finally interacts with a nucleus.

What happens when it interacts is then decided by the use of further random numbers to select the type of nuclide it interacts with and the interaction (elastic scattering, inelastic scattering, capture or fission) that takes place, the choices being weighted by the various macroscopic cross-sections. In the case of elastic or inelastic scattering the neutron continues on its way in a different direction and with a different energy, again chosen at random taking account of any possible anisotropy of elastic scattering and the probabilities of exciting different nuclear energy levels in the case of inelastic scattering. In the case of capture the neutron disappears. In the case of fission it is replaced by one or more new neutrons, their number and energies chosen at random in accordance with data on v and x for the fissile nuclide concerned. These new neutrons are then tracked as before and the tracking continues until all daughter neutrons are captured or have travelled out of the reactor.

The same procedure is followed for many neutrons and accounts are kept of the numbers and locations of the different interactions. For example reaction rates can be estimated by keeping account of the number of reactions taking place in an element of volume, and the power density can be deduced from the rate at which fissions take place. Neutron flux can be estimated by keeping account of the total path-length traversed by neutrons in the volume element, and neutron current by noting the number of neutrons that pass through an element of area. An estimate of one of these quantities is given by limN^TO (Nq/N), where Nq is the number of records of event q and N is the total number of neutron histories that have been calculated. Additional neutrons are tracked until a statistical test shows that the variance of the estimates of q is small enough to give confidence that

the limit has been approached sufficiently closely. The most intricate part of the calculation is the determination of which regions a neutron passes through and where it crosses the inter-region boundaries.

The procedure for a critical reactor is a little more complicated than that for one driven by a neutron source because as well as the fluxes and reaction rates the eigenvalue ke has to be estimated. It can be done as follows. A number N1 of neutrons, sometimes called a “batch”, are tracked up to the point at which they are captured, leak out of the reactor, or cause fissions. If a total of N2 new neutrons are born in these fissions, N2/N1 = ke1 is a first estimate of ke. All values of v are then divided by this estimate and the N2 new neutrons from the fissions caused by the first batch, which become the second batch, are tracked until they in turn cause fissions that give rise to N3 further neutrons. (N3/N2)ke1 = ke2 is then a new estimate of ke, the v are modified again, a third batch of neutrons are tracked, and so on. In due course the estimates cluster around a limit, which is the required value of ke, and the process continues until the variance is small enough. In making the final estimate of ke it is normal to reject the values produced by the initial batches of neutrons as they converge towards the final value. The process is shown diagrammatically in Figure 1.4. The neutron fluxes generated in the course of this process are of course meaningless because they relate to a reactor with artificial values of v. It is necessary to iterate by adjusting its composition or dimensions until ke = 1.

To obtain reliable estimates of ke or the critical composition it is necessary to track large numbers of neutrons. For a large power reactor a batch would consist of 104 or 105 neutrons and several hun­dred batches would be calculated with possibly the first 50 or so being rejected as the estimates converged. In all tens of millions of neutron histories are likely to be followed to give an estimate of ke with a stand­ard deviation of 0.0001 or better. No doubt the numbers of neutron histories will increase and the ke estimates will improve as computing facilities become more powerful. (The numbers can be put in perspect­ive by remembering that a 2500 MW (heat) reactor is producing some 2 x 1020 fission neutrons per second.)

The variance of the estimate of a quantity such as the flux at a certain point in the reactor depends on the number of neutron histor­ies that contribute to it. This implies that more histories have to be calculated for a reliable estimate of the flux in a low-flux region than for the flux at the centre of the core. Since the variance is roughly inversely proportional to the square root of the number of histories this implies that if it is necessary to calculate N histories to attain a satisfactory estimate of the flux at the core centre, 100 N histories will be required for an equally reliable estimate of the flux in a region, such as a breeder or a reflector, where it is a tenth of that at the centre. This consideration would, in principle, make it impossible to use a Monte Carlo method to calculate the performance of a reactor shield, which may be required to attenuate the flux by a factor of 1012 or more, but the difficulty can be avoided by “splitting” neutrons.

For example if a neutron crosses the outer boundary of the reactor core and enters the shield it can be replaced by two neutrons each of which has half the “weight” of the original neutron and is then allowed to have an independent history. The process can be repeated at other boundaries in the outer parts of the shield so that there are statistically

significant numbers of neutrons throughout the shield but the atten­uation is accounted for by their diminished weights. The process is shown diagrammatically in Figure 1.5.