The Critical Mass

One of the questions surrounding the Comprehensive Test Ban Treaty is what can an opponent gain by cheating. It may be possible (however unlikely) to set off small nuclear explosions that are below the detection threshold of about 1 kiloton for the International Monitoring System. Large underground cavities are used to decouple the blast from the surrounding rock and suppress the seismic signature of the explosion. Important technical questions associated with this scenario are what can one learn from such tests and can that knowledge be extrapolated to higher yields? The answers to these questions depend on the quality of the physics models used to describe the process. In this computational laboratory we will explore a simple model of neutron diffusion. Nuclear fission is caused by absorption of a neutron by a fissile nucleus like 235-U and 239-Pu so understanding how the cloud of neutrons spreads through the fissile material is essential to understanding the physical process. The validity of these models is essential for the policy questions raised above. The goal of this laboratory is to find the critical mass or size of a one-dimensional 'pipe' of 235-U that will explode.

1.0 Some Necessary Commands

In this section several Mathematica commands that will be useful in your calculation are demonstrated.  You should execute the commands in the cells below to investigate the new functions.

1.1 Explicit Methods for Solving Partial Differential Equations.

The diffusion equation in one dimension and for many cases can be written in the following way

((∂n(x,t))/(∂ t))=D((((∂ )^2)n(x,t))/(∂ (x^2)))

where n(x,t) is the neutron density as a function of position and time and D is the coefficient of self-diffusion. In this section we demonstrate a code that uses an explicit method to solve the diffusion equation.

(* Define diffusion parameters. *)
Dn=0.001;  (* self diffusion coefficient in (m^2)/shake *)
Ln = 0.10;  (* size of the region in meters. *)

(* parameters for the algorithm. *)
tmax = 10.0;  (* maximum time in shakes. *)
Nxsteps = 20;  (* steps in x. *)
Ntsteps = 1000;  (* steps in time. *)
dx = (Ln/Nxsteps);  (* stepsize in x (m). *)
dt = (tmax/Ntsteps);  (* stepsize in time. *)

(* set up the distribution of particles at t=0. *)
n0=Table[{x,0,0},{x,0,Ln,dx}];
n0[[(Nxsteps/2)+1,3]] = (1/dx) ;

(* initialize the main array. *)
particle = Table[0.0,{i,1,Nxsteps},{n,1,Ntsteps}];
(* put in the initial conditions for t=0. *)
Do[particle[[i,1]]=n0[[i,3]],{i,1,Nxsteps}];
(* The boundary condition at x=0. *)
Do[particle[[1,n]]=10.0,{n,2,Ntsteps}];
(* The boundary condition at x=L. *)
Do[particle[[Nxsteps,n]]=0.6,{n,2,Ntsteps}];

(* constants for the recursion relation. *)
A0 = 1 - ((2*dt*Dn)/(dx^2));
B0 = ((dt*Dn)/(dx^2));

(* main loop. outer loop over time and inner loop over position. *)
Do[
Do[particle[[i,n]] = A0*particle[[i,n-1]] + B0*particle[[i+1,n-1]] + B0*particle[[i-1,n-1]],
{i,2,Nxsteps-1}
],
{n,2,Ntsteps}
]

1.2 2-D Plotting.

It is sometimes useful to plot the particle density as a log-log plot since things happen over a very large dynamic range.

(* Pick out the x-value in the middle of the range. *)
xcounter=IntegerPart[Nxsteps/2];
xvalue=dx*xcounter;
t1=Table[{dt*(n-1),particle[[xcounter,n]]},{n,2,Ntsteps}];
ListLogLogPlot[t1,PlotRange->Automatic,Frame->True,FrameLabel->{"t (s)","n ((m^(-1)))",StringForm["Diffusion for x=`` m",xvalue],""},Joined->True,
BaseStyle->Large,PlotStyle->Thickness[0.005],LabelStyle->Directive[Larger]]

1.3 3-D Plotting.

Since our particle density n(x,t) is three dimensional, it helps to plot it like that. Take care because these plots can take a long time to complete.

(* Create a 3D table of position, time, and particle density. *)
t2=Table[{dx*(i-1),dt*(n-1),particle[[i,n]]},{i,1,Nxsteps,1},{n,1,Ntsteps,1}];
(* The previous Table command puts brackets around each set of position points associated with a specific value of time. We need to remove those brackets to make the plot we want. *)
t2b=Flatten[t2,1];
ListPlot3D[t2b,AxesLabel->{"x (m)","t (s)"," n ((m^(-1)))"},PlotRange->{0,20}]

2.0 Will It Go Boom?

We want to explore the conditions under which a one-dimensional pipe of 235-U will explode. What this means in our model is that the neutron density will rise rapidly with time. The process starts when two pieces of 235-U come together and a neutron fuse is set off. In our model we have two, identical, long uniform pipes of fissile material that are brought together rapidly. A neutron trigger at the center causes the fission reaction to start. 

a. The test code in Part 1 included a spike in the particle density at the center of the pipe at t=0. Use this same configuration for your code.

b. Modify the code in Part 1 to include a source term (C≠0) in the diffusion equation.

c. The parameters for 235-U are the following. You can choose the length L of the pipe to be whatever you want.

D=(10^5) (m^2)/s
C=(10^8) (s^(-1))

The boundary conditions are that neutrons that reach the ends of the pipe leave forever so require that the neutron density n(x,t) go to zero at the ends of the pipe.

d. Under what conditions will your 'bomb' explode? Why?