Student-level numerical simulation of conditions inside an exploding fission-bomb core

This paper describes a freely-available spreadsheet that has been developed to simulate the conditions of reaction rate, core acceleration and velocity, energy generation, and pressure within a detonating fission-bomb core. When applied to a model of the Hiroshima Little Boy bomb, the spreadsheet predicts a yield of 12.7 kilotons, a figure in reasonable agreement with published values.


INTRODUCTION
The discipline of computational physics is now regarded as possessing importance equal to the traditional areas of experimental and theoretical studies. As such, it is important that our students be introduced early in their careers to the power of modern desktop computational tools that can be used to model physically interesting systems. Computational models can facilitate study of physical systems where the theory cannot be solved purely analytically and/or that are not easily realizable experimentally. Sometimes a single graph can serve to dramatically make a point about the behavior of a system that is "latent" in the mathematics but whose magnitude is not immediately apparent. This paper describes and makes freely available an Excel spreadsheet to carry out a student-level simulation of a system that exemplifies all of these points: the conditions of energy release, pressure, fission rate, and expansion inside the core of a detonating fission bomb. The physics of nuclear weapons has been and will remain a fascinating and timely subject. These devices are forbidding and mysterious to students and the general public alike; by helping our students to understand some of the details of their functioning we can equip them to play constructive roles in furthering public understanding of them.
The basic physics of criticality conditions in both "bare" (untamped) and tamped fission cores is described in two papers previously published by this author [1,2]. Because the time-dependence of conditions within a fissioning core is highly non-linear, approximate analytic expressionss had to be developed in those papers in order to arrive at estimates of bomb yields and efficiencies. The purpose of this paper is to utilize the theoretical foundations established in those papers to build an approximate but easy-to-use numerical simulation of a core whose mass, nuclear properties, and tamper properties are set by the user. The resulting spreadsheet is made freely downloadable to interested readers. The structure of this paper is as follows. In Section 2 I lay out the theoretical background; while this is adopted directly from References [1] and [2] it is presented here for sake of a self-contained discussion. The programming of the simulation itself is described in Section 3, and in Section 4 I present the results of a simulation of the Hiroshima Little Boy uranium bomb. A brief summary is presented in Section 5.

FISSION-CORE PHYSICS
The essential physical quantity within a fissioning bomb core is the number density of neutrons, N. As described in References [1] and [2], diffusion theory leads to an approximate expression for the space and time-dependence of N within a spherically symmetric core of the form In this expression, N o is the initial neutron number density at the center of the core, and is set by whatever initiating device starts the chain reaction. is the average time that a neutron will travel before causing a fission, neut core where v neut is the average neutron speed and core fiss core, where n is the nuclear number density and f  the fission cross-section of the fissile core material. d core in Eq.1 has units of length and can be thought of as fundamentally setting the size of the critical radius; it is given by where  is the number of neutrons emitted per fission (so-called secondary neutrons) and  trans core is the transport mean free path for neutrons,  t is the so-called transport cross-section. If non-fission neutron capture can be ignored (which should be the case for any sensibly pure fissile material), this is given by the sum of the fission and elastic-scattering crosssections: The parameter  in Eqs.1 and 4 arises in a separation of variables in solving the diffusion equation for the neutron density, and is itself time-dependent as it depends on the core radius as described in what follows. First consider the case of an untamped core, also known as a bare core. As described in References [1] and [2], application of appropriate boundary conditions to Eq.1 to ensure that there is no way for neutrons which have passed out through the surface of the core to be reflected back inside it from the external world leads to the following constraint for the core radius R C : This is the criticality condition for a bare core. Satisfaction of this constraint depends on through its appearance in d core . If = 0, then the neutron number density is neither growing nor diminishing in time, a condition known as threshold criticality. (In an implosion weapon, the moment when this state first occurs during compression of a subcritical core into a supercritical mass is known as first criticality. The present work is not inteneded to model implosion scenarios, which are much more complex on account of the inward motion of the core at the moment when fissions are initiated.) In this case d core is completely determined by the fissility parameters, and the value of R C which satisfies Eq.7 is the bare threshold critical radius, thresh bare R ; for pure U-235 this is about 8.4 cm, equivalent to about 46 kg.
As explained in Reference [1], if one starts with a core of specified radius R C > thresh bare R , then Eq.7 can be solved for  (through d core ); one will find that > 0, which means that the neutron density is growing exponentially in time, a condition known as "supercriticality." As the core rapidly (within a microsecond) expands due to the extreme rate of release of energy by fissions, will decline as a function of time until it reaches zero, at which point the chain reaction will rapidly shut down; this situation is known as "second criticality" and effectively marks the end of the detonation phase.
In Reference [2], it is shown that if the core is surrounded by a snugly-fitting tamper of non-fissile and non-neutron-absorbing material, diffusion theory leads to the following expressions for the neutron density within the tamper material: where  is as above and  is again the mean travel time for neutrons between fissions in the core, A and B are constants of integration to be determined by boundary conditions, and d tamp is given by where tamp trans  is the mean free path for neutron transport in the tamper material, analogous to Eq.5 except that the tamper material is presumed to have no fission crosssection: If the core and tamper have outer radii R core and R tamp , then demanding the continuity of neutron density and flux at the core/tamper interface and again requiring that that no neutrons which escape from the tamper to the external world can be reflected back, one finds that the criticality conditions emerge as and Eq.11 corresponds to tamped threshold ( = 0) criticality. Once values for the d's and 's are given the only unknown is R core thresh , the core radius for tamped threshold criticality. In using Eqs.12 and 13 the idea is that one again specifies the mass and hence radius of the core, R core (> thresh core R ), and solves for . This would presumably be the value of  when fissions are initiated at t = 0; as the core expands  will subsequently decline until second criticality is reached. A tamper serves to increase the efficiency (and hence the yield) of a weapon through two effects. First, by briefly retarding the expansion of the core, the tamper causes  to remain greater for a longer time than it would have otherwise; this is beneficial as the rate of fissions -and hence the energy release -depends exponentially on . Second, the tamper serves to reflect neutrons back into the core, effectively decreasing the loss of fission-causing neutrons from within the core to the outside world. In modern weapons -engineering parlance a tamper is known as a reflector; I retain the historical terminology. The retardation effect is difficult to model analytically at this level and so is not accounted for in Eqs.12 and 13; they do, however, include the reflection effect in the boundary conditions used to establish them. The retardation effect is treated approximately in the simulation as described following Eq.17 below.
We can now consider the time-dependence of various quantities in order to begin formulating a simulation. At a moment when the core has volume V core , the rate of fissions R(t) (fission/sec) is given by If each fission releases energy E fiss (typically ~ 180 MeV), then the rate of energy release within the core is The total energy liberated to a given time can be tracked by numerically integrating Eq.15; this determines the pressure within the core as a function of time. This follows from the thermodynamic pressure -energy density relationship The choice of the parameter  depends on whether gas pressure ( = 2/3) or radiation pressure ( = 1/3) is dominant; the latter dominates for per-particle energies greater than about 2 keV and will presumably be the case for the later, more energetic stages of the reaction. I use the core volume in Eq.16 on the rationale that the fission products which cause the gas/radiation pressure will likely largely remain within the core. Following Reference [1], I model the core as an expanding sphere of radius r(t) with all parts of the sphere moving at speed v(t), driven by the energy release from fissions. Do not confuse this velocity with the average neutron speed, which does not directly come into this part of the development (it does enter implicitly, however, through ). Invoking the work-energy theorem in its thermodynamic formulation, W = P(t) dV, I equate the work done by the gas (or radiation) pressure in changing the core volume by dV over time dt to the change in the core's kinetic energy over that time: For simplicity in developing the simulation, I treat the tamper as remaining of constant density. Now, it is desirable to make some effort to account for the retarding effect of the tamper on the core. To do this, I treat the dK/dt term in Eq.17 as involving the velocity of the core expansion but with the mass involved being that of the core plus that of the tamper. The dV/dt term is taken to apply to the core. With r as the radius and v the expansion velocity of the core, we have Numerically integrating this result will give v(t) = dr/dt, which can be integrated to give the core radius as a function of time.

3.THE SIMULATION
I have developed an Excel spreadsheet to carry out the calculations described above. This is freely available at http://othello.alma.edu/~reed/FissionCore.xls. In this Section I describe the general layout of this spreadsheet; some results are described in Section 4 below. This spreadsheet consists of three interlinked sheets. On the first, the user inputs fundamental data such as core and tamper material densities, atomic weights, cross-sections, the secondary neutron number, the average secondary-neutron energy, values for E f and , the desired core mass, the outer radius of the tamper, and the number of "initial neutrons" in the core at t = 0. These are entered in convenient units such as g/cm 3 , barns, and MeV; the spreadsheet subsequently carries out all calculations in MKS units. The Excel "Goal Seek" function is then run three times, to establish values for 1) the bare threshold critical radius, 2) the tamped threshold critical radius, and 3) the value of  corresponding to the chosen core mass. The masses in 1) and 2) are computed for reference and for the fact that they are needed for some calculations involving the expansion of the core as described below. The chosen core mass should exceed that corresponding to thresh tamp R .
A significant complexity in carrying out this simulation is that one apparently needs to solve Eq.12 for the value of  corresponding to each time-stepped core radius between first and second criticality: the fission rate, energy generation rate, and pressure all depend on  as a function of time. I have found, however, that  is usually quite linear as a function of core radius. This behavior greatly simplifies the actual time-dependent simulation. Sheet 2 of the spreadsheet allows one to establish parameters for this linear behavior for the values of the various parameters that the user inputs on Sheet 1. Here, the user solves (again using the Goal Seek function) for the value of  for 25 values of the radius. These start at the initial core radius and proceed to 1.25 times the value of the second-criticality radius for a bare core of the mass chosen by the user on Sheet 1; this range appears to be suitable to establish the behavior of . The rationale for this arrangement is as follows. As shown in Reference [2], if the chosen core mass is equal to C bare threshold critical masses, criticality will hold over a range of radii given by   The presence of a tamper means that the core will expand somewhat beyond r before second criticality is reached, but Eq.19 sets the essential length scale of the expansion. For convenience, Sheet 2 utilizes a "normalized" radius defined as where C is now defined as the number of tamped threshold critical masses. r norm = 1 corresponds to the second criticality radius one would compute from Eq.19 if it applied as well to a tamped core. Sheet 2 tracks the changing mass density, nuclear number densities, and mean-free-paths within the core as a function of r. By running the Goal Seek function on each of the 25 radii, the user adjusts  in each case to render Eq.12 equal to zero. The behavior of (r) is then displayed in an automatically-generated graph. On a separate line with fixed to a value very near zero (10 -10 is built-in), the user adjusts the radius to once again render Eq.12 equal to zero, thus establishing the radius of second criticality for his or her parameters. The slope and intercept of a linear (r) fit are then automatically computed in preparation for the next step. The actual time-dependent simulation occurs on Sheet 3. The simulation is set up to involve 500 timesteps, one per row. The initial core radius is transferred from Sheet 1 for t = 0. Because much of the energy release in a nuclear weapon occurs during the last few generation of fissions before second criticality, this Sheet allows the user to set up two different timescales: an "initial" one (dtinit) intended for use in the first few rows of the Sheet when a larger timestep can be tolerated without much loss of accuracy, and a later one (dtlate), to be chosen considerably smaller and used for the majority of the rows. In this way a user can optimize the 500 rows to both capture sufficient accuracy in the last few fission generations and arrange that (r) is just approaching zero at the last steps of the process. Typical choices for dtinit and dtlate might be a few tenths of a microsecond and a few tenths of a nanosecond, respectively. At each radius, the Sheet computes the value of (r) from the linear approximation of Sheet 2, the core volume, mass density, nuclear number densities and mean free paths within the core, , rates of fission and energy generation, pressure, and total energy liberated to that time. The acceleration of the core is computed from Eq.18, and the core velocity and radius are updated depending upon the timestep in play; the new radius is transferred to the subsequent row to seed the next step. The user is automatically presented with graphs of (r), the fission rate, pressure, and total energy liberated (in kilotons equivalent) as functions of time.

A SIMULATION OF THE HIROSHIMA LITTLE BOY BOMB
As described in References [2] and particularly [3], the Hiroshima Little Boy core comprised about 64 kg of enriched U-235 in a cylindrical configuration surrounded by a cylindrical tungsten-carbide tamper of diameter and length 13 inches, mass approximately 310 kg, and density 14.8 g /cm 3 . Values for the various core and tamper parameters are given in Table 1; these are adopted from Reference [2]. Assuming these values and taking the core to be spherical (radius 9.35 cm at a density of 18.71 g/cm 3 ; this figure is 235/238 times the density of natural uranium, 18.95 g/cm 3 ) and surrounded by a spherical tungsten-carbide tamper of outer radius 18 cm (mass 311 kg), Sheet 1 of the author's spreadsheet indicates that the tamped threshold critical mass of U-235 in this configuration is 18.4 kg, about a 60% reduction from the bare threshold critical mass of 45.9 kg. Figure 1 shows that the run of (r) for this situation is quite linear out to the computed second criticality radius of 12.04 cm. Upon Table 1. Parameters for U-235 core/tungsten-carbide tamper model.   within an interval of about 0.1 s. The pressure peaks at close to 5 × 10 15 Pa, or about 50 billion atmospheres, equivalent to about one-fifth of that at the center of the Sun, and the fission rate peaks at about 4 × 10 31 per second. The core acceleration peaks at about 1.5 × 10 12 m/s 2 at t -0.9 s, and second criticality occurs at t -1.05 s, at which time the core expansion velocity is about 270 km/s. These graphs dramatically illustrate what Robert Serber wrote in The Los Alamos Primer : "Since only the last few generations will release enough energy to produce much expansion, it is just possible for the reaction to occur to an interesting extent before it is stopped by the spreading of the active material" [4]. The predicted yield of Little Boy from the present model is 12.7 kt. This result is in surprisingly good agreement with the estimated -12 kt yield published by Penney, et al. [5]. At a fission yield of 17.59 kt per kg of pure U-235, this represents an efficiency of only about 1.1% for the 64 kg core. Some of this agreement must be fortuitous, however, in view of the approximations incorporated in the present model. That the yield estimate needs to be taken with some skepticism is demonstrated by the fact that increasing the initial number of neutrons to 10 increases the yield to 18.7 kt. However, this change does not much affect the timescale or the peak pressure and fission rates. Users who download FissionCore.xls will find this example pre-loaded.

SUMMARY
This paper describes the development of a spreadsheet for simulating the conditions within a detonating fission-bomb core. The simulation is straightforward enough to be used with students, and for a simulation of