Project 03: What volumn of an impure $^{235}\text{U}$ that yields a steady-state behaviour in the 3rd generation nuclear chain reaction?¶
Name: Nam Bui
Introduction:¶
What volume of an impure $^{235}\text{U}$ that yields a steady-state behaviour in the 3rd generation nuclear chain reaction under these conditions:
- 𝑁0=250 initial thermal neutrons randomly distributed inside the cubic atom.
- Neutrons will randomly scatter with a mean free path of 𝜆=0.15 m.
- If neutrons ended up inside the atom, they would have a chance to undergo nuclear fission, creating two new neutrons.
- The probability of fission is lower near one vertex of the atom (origin) and higher the further from the origin. $$p = \frac{\sqrt{x_0^2 + y_0^2 + z_0^2}}{a}$$, where a is the side length of the cube. The reason for this definition is that it is easier for me to formulate and simulate if choosing the distance from the origin.
- The side length (a) will be varied, thus the volume in order to find out which volume would yield a stable behaviour k = 1 in the 3rd generation.
- The distance travelled by each particle comes from the probability, 𝑝(𝐿)∝(1/λ)𝑒−𝐿/λ which I will use d = np.random.exponential(mean_free_path, N0) in my code.
- The steady-state behaviour of the system is indicated by the k factor, which is the ratio of number of particles between two generations when a new generation is created.
- 𝑘<1 : Subcritical (the reaction dies out); 𝑘=1 : Critical (the reaction is self-sustaining); or 𝑘>1: Supercritical (the reaction grows exponentially).
- There will be 250 trials to generate the population of each generation which is then used to determine the k factor. The mean k factor of each generation will be calculated and used for the report.
The Result:¶
import numpy as np
import matplotlib.pyplot as plt
V_list = []
k23_mean_list = []
for a in np.arange(0.1, 0.5, 0.04):
V = a**3
V_list.append(V)
# Parameters
N0 = 250 # Initial number of neutrons
mean_free_path = 0.15 # Mean free path (m)
m_trials = 250 # Number of trials
k23_values = [] # Store k23 for each trial
for trial in range(m_trials):
neutron_fiss_1 = 0
positions_1 = []
d = np.random.exponential(mean_free_path, N0)
for dist in d:
x0, y0, z0 = np.random.uniform(0, a, size=3)
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_1 += 1
positions_1.append((x0, y0, z0))
N1 = neutron_fiss_1 * 2
neutron_fiss_2 = 0
positions_2 = []
d2 = np.random.exponential(mean_free_path, N1)
for dist, (x0, y0, z0) in zip(d2, positions_1 * 2):
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_2 += 1
positions_2.append((x0, y0, z0))
N2 = neutron_fiss_2 * 2
neutron_fiss_3 = 0
d3 = np.random.exponential(mean_free_path, N2)
for dist, (x0, y0, z0) in zip(d3, positions_2 * 2):
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_3 += 1
N3 = neutron_fiss_3 * 2
k23 = N3 / N2
k23_values.append(k23)
k23_mean = np.mean(k23_values)
k23_mean_list.append(k23_mean)
# Find the volume that yields k = 1
stable = [(i, val) for i, val in enumerate(k23_mean_list) if abs(val - 1) <= 0.02]
stable_volume = V_list[stable[0][0]] if stable else None
if stable_volume:
print(f"The volume that yields k = 1 is: {stable_volume:.5g} m^3")
# Plot
plt.figure(figsize=(8, 6))
plt.plot(V_list, k23_mean_list, marker='o', linestyle='-', label='Mean k23 vs Volume')
if stable_volume:
plt.axvline(x=stable_volume, color='r', linestyle='--', label=f'Stable Volume: {stable_volume:.3g} m^3')
plt.xlabel('Volume (m^3)')
plt.ylabel('Mean k23')
plt.title('Relationship Between Atom Volume and k23')
plt.grid(True)
plt.legend()
plt.show()
The volume that yields k = 1 is: 0.039304 m^3
Conclusion:¶
The simulation found out that the impure cubic $^{235}\text{U}$ atom with volume: 0.039304m^3 will yield a stable k factor in the third generation of neutrons.
Indeed, this is restricted to a cubic-shaped atom only and the mean free path is restricted to 0.15m. Another restriction of this simulation is the initial population of neutrons is not quite big enough (250 particles), and I could also increase the trials too. Otherwise, the simulation helped answer my research question effectively.
To leverage my simulation, I will try to shorten the running time by minimizing "for" loop in order to look at a larger neutron initial population as well as trials.
References¶
(Each reference in this section should include a citation associated with that reference in the body of your project. You may find it helpful to look at a typical Wikipedia page as an example of how to do this. Make sure you include references for source of datas, literature cited, and any python packages used beyond our standard ones. Additionally, you must also cite the sources for any code that you found on the internet or from peers.)
Appendix 1: Code validation¶
A1.1: Validations of Position Distribution and Distance Travelled From Two Values of "a"¶
In this Appendix, I will examine one volume with a stable 3rd generation k (a=0.34) and one with unstable (a=0.46) to see if neutrons' random movement and distribution were obeyed. In addition, I'm going to examine if the distance travel of neutrons actually followed the mean free path and obeyed the exponential distribution of 𝑒^−𝐿/λ. In the code below, I will generate a histogram of the distance travelled by neutrons for the first trial only and extract their mean value. Also, I will plot the positions of neutrons after their first random movement for examination.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
V_list = []
k23_mean_list = []
labels = []
for a in [0.34, 0.46]: # Loop through both cube sizes
V = a**3
V_list.append(V)
# Add labels for stability
if a == 0.34:
labels.append("Stable k (a=0.34)")
else:
labels.append("Unstable k (a=0.46)")
# Parameters
N0 = 250 # Initial number of neutrons
mean_free_path = 0.15 # Mean free path (m)
m_trials = 250 # Number of trials
k23_values = [] # Store k23 for each trial
for trial in range(m_trials):
# First generation
neutron_fiss_1 = 0
positions_1 = [] # List to store positions where first-generation neutrons stop
# Generate travel distances for all neutrons in N0
d = np.random.exponential(mean_free_path, N0)
x0, y0, z0 = np.random.uniform(0, a, size=(3, N0))
for i, dist in enumerate(d):
# Generate random direction
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
# Final position after traveling distance d
x1 = x0[i] + dist * dx
y1 = y0[i] + dist * dy
z1 = z0[i] + dist * dz
# Probability of fission:
p = np.sqrt(x1**2 + y1**2 + z1**2) / a
# Check if neutrons undergo fission:
if 0 <= x1 <= a and 0 <= y1 <= a and 0 <= z1 <= a:
if np.random.rand() < p:
neutron_fiss_1 += 1
positions_1.append((x1, y1, z1)) # Save stopping position
# Validation 1: Visualizing travel distances and positions (for the first trial only)
if trial == 0:
# Plot histogram of travel distances
plt.figure(figsize=(8, 6))
plt.hist(d, bins=20, alpha=0.7, density=True, label=f"Travel distances for a={a}")
# Add mean distance
mean_distance = np.mean(d)
plt.axvline(mean_distance, color='r', linestyle='--', label=f"Mean distance: {mean_distance:.3f} m")
# Finalize plot
plt.xlabel("Distance (m)")
plt.ylabel("Probability Density")
plt.title(f"Distribution of Neutron Travel Distances (a={a})")
plt.legend()
plt.show()
# 3D scatter plot of neutron positions
positions_1_array = np.array(positions_1)
if len(positions_1_array) > 0: # Ensure there are positions to plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(
positions_1_array[:, 0],
positions_1_array[:, 1],
positions_1_array[:, 2],
alpha=0.6,
label=f"Neutron positions (a={a})",
c=positions_1_array[:, 2], # Color by z-coordinate for depth effect
cmap="viridis"
)
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
ax.set_zlabel("z (m)")
ax.set_title(f"Neutron Positions in the Cube (a={a})")
plt.legend()
plt.show()
# Calculate N1 and k01
N1 = neutron_fiss_1 * 2
# Second generation
neutron_fiss_2 = 0
positions_2 = [] # List to store positions where second-generation neutrons stop
d2 = np.random.exponential(mean_free_path, N1)
for dist, (x0, y0, z0) in zip(d2, positions_1 * 2): # Each position in positions_1 produces 2 new neutrons
# Generate random direction
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
# Final position after traveling distance d
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
# Probability of fission:
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
# Check if neutrons undergo fission:
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_2 += 1
positions_2.append((x0, y0, z0)) # Save stopping position
# Calculate N2 and k12
N2 = neutron_fiss_2 * 2
# Third generation
neutron_fiss_3 = 0
d3 = np.random.exponential(mean_free_path, N2)
for dist, (x0, y0, z0) in zip(d3, positions_2 * 2): # Each position in positions_2 produces 2 new neutrons
# Generate random direction
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
# Final position after traveling distance d
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
# Probability of fission:
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
# Check if neutrons undergo fission:
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_3 += 1
# Calculate N3 and k23
N3 = neutron_fiss_3 * 2
k23 = N3 / N2
k23_values.append(k23)
# Calculate averages and standard deviations for k23
k23_mean = np.mean(k23_values)
k23_mean_list.append(k23_mean)
Indeed, from the plots, I could see that both "a" values yielded a uniform random distribution of neutrons inside the cubic atom even after the neutrons randomly travelled once (first generation). Also, the mean distance travel (0.146m) of neutrons obeyed the mean free path (0.15m), as both values are close.
A1.2: Fission Probability Validation:¶
In this Appendix, I will look at the fissioned neutrons in the 2nd generation for the first trial only. Therefrom, I could see the pattern if the fission really followed my definition that the probability is less near (0,0,0) and more the further away.
import numpy as np
import matplotlib.pyplot as plt
V_list = []
k23_mean_list = []
for a in [0.34,0.46]:
V = a**3
V_list.append(V)
# Parameters
N0 = 250 # Initial number of neutrons
mean_free_path = 0.15 # Mean free path (m)
m_trials = 250 # Number of trials
k23_values = [] # Store k23 for each trial
for trial in range(m_trials):
# First generation
neutron_fiss_1 = 0
positions_1 = [] # List to store positions where first-generation neutrons stop
# Generate travel distances for all neutrons in N0
d = np.random.exponential(mean_free_path, N0)
for dist in d:
# Generate initial position
x0, y0, z0 = np.random.uniform(0, a, size=3)
# Generate random direction
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
# Final position after traveling distance d
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
# Probability of fission:
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
# Check if neutrons undergo fission:
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_1 += 1
positions_1.append((x0, y0, z0)) # Save stopping position
# Calculate N1 and k01
N1 = neutron_fiss_1 * 2
# Second generation
neutron_fiss_2 = 0
positions_2 = [] # List to store positions where second-generation neutrons stop
fission_flags_2 = [] # To store fission status of second-generation neutrons
d2 = np.random.exponential(mean_free_path, N1)
for dist, (x0, y0, z0) in zip(d2, positions_1 * 2): # Each position in positions_1 produces 2 new neutrons
# Generate random direction
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
# Final position after traveling distance d
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
# Probability of fission:
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
# Check if neutrons undergo fission:
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_2 += 1
fission_flags_2.append(True) # Mark as undergoing fission
else:
fission_flags_2.append(False) # Mark as not undergoing fission
positions_2.append((x0, y0, z0)) # Save stopping position
# Plot 3D positions for the second generation (only for the first trial)
if trial == 0:
positions_2_array = np.array(positions_2)
fission_positions = positions_2_array[np.array(fission_flags_2)]
no_fission_positions = positions_2_array[~np.array(fission_flags_2)]
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# Plot neutrons that undergo fission
ax.scatter(
fission_positions[:, 0], fission_positions[:, 1], fission_positions[:, 2],
c='red', alpha=0.6, label="Neutrons undergoing fission"
)
# Plot neutrons that don't undergo fission
ax.scatter(
no_fission_positions[:, 0], no_fission_positions[:, 1], no_fission_positions[:, 2],
c='blue', alpha=0.6, label="Neutrons not undergoing fission"
)
# Set labels and title
ax.set_xlabel("x (m)")
ax.set_ylabel("y (m)")
ax.set_zlabel("z (m)")
ax.set_title(f"Second-Generation Neutron Positions (a={a:.2f})")
ax.legend()
plt.show()
# Calculate N2 and k12
N2 = neutron_fiss_2 * 2
# Third generation
neutron_fiss_3 = 0
d3 = np.random.exponential(mean_free_path, N2)
for dist, (x0, y0, z0) in zip(d3, positions_2 * 2): # Each position in positions_2 produces 2 new neutrons
# Generate random direction
phi = np.random.uniform(0, 2 * np.pi)
costheta = np.random.uniform(-1, 1)
theta = np.arccos(costheta)
dx = np.sin(theta) * np.cos(phi)
dy = np.sin(theta) * np.sin(phi)
dz = np.cos(theta)
# Final position after traveling distance d
x0 += dist * dx
y0 += dist * dy
z0 += dist * dz
# Probability of fission:
p = np.sqrt(x0**2 + y0**2 + z0**2) / a
# Check if neutrons undergo fission:
if 0 <= x0 <= a and 0 <= y0 <= a and 0 <= z0 <= a:
if np.random.rand() < p:
neutron_fiss_3 += 1
# Calculate N3 and k23
N3 = neutron_fiss_3 * 2
k23 = N3 / N2
k23_values.append(k23)
# Calculate averages and standard deviations for k23
k23_mean = np.mean(k23_values)
k23_mean_list.append(k23_mean)
From the plots, I could see that the not-undergo fission (blue dots) are populated more near x,y,z = 0 than the other places. Thus, the simulation worked as I intended.
A1.3:¶
Appendix 2: Reflection questions¶
Reflection 1: Coding Approaches (A)¶
(How well did you apply and extend your coding knowledge in this project? Consider steps you took to make the code more efficient, readable and/or concise. Discuss any new-to-you coding techniques, functions or python packages that you learned how to use. Reflect on any unforeseen coding challenges you faced in completing this project.)
In my opinion, I effectively applied coding techniques taught in class, as well as new techniques I learned from online sources and generative AI. Typical Joss's techniques that I applied are np.random.exponential and np.random.uniform for the random distance travel that saturates around the mean free path and the uniform random motion of each particle; the Monte Carlo's technique for deciding if a particle will undergo fission using np.random.rand() < p. Besides, I learned new list comprehension from generative AI for detecting the volume that yields a stable k factor. It is new to me since AI uses a tuple of two variables in the list. One challenge regarding my coding process is to detect the volume that yields a stable k factor since I need to choose a step size small enough so that the k factor is within 0.02 units from 1 (stable value), while the running time for my code is quite long.
Reflection 2: Coding Approaches (B)¶
(Highlight an aspect of your code that you feel you did particularily well. Discuss an aspect of your code that would benefit the most from further effort.)
The part of the code that I found I did well is how I found and marked in the plot the neutron's positions which undergo fission. I first flag each neutron's position as "True" or "False", then I divided them into two lists (fission and not fission) using the flags. I also learned from gen AI to use ~ (tilde) instead of using list comprehension "if not..." for the list dividing task. The aspect of the code that I think would benefit most if I put more effort is to improve the overuse of the "for" loops which makes my code run very slow.
Reflection 3: Simulation phyiscs and investigation (A)¶
(How well did you apply and extend your physical modelling and scientific investigation skills in this project? Consider the phase space you chose to explore and how throroughly you explored it. Consider how you translated physics into code and if appropriate any new physics you learned or developed a more thorough understanding of.)
For this project, I think I did quite well in incorporating the essential physical models: mean free path, my made-up fission probability, and uniform random motion using vectors. For the mean free path, it is actually quite simple to include since it is available in numpy library. For the random uniform motion, Joss's given codes provide the unit vector in 3 directions x,y, and z translated from the spherical coordinate (phi angle from x axis and theta as angle from z axis). Then I only needed to randomize the unit vectors and multiply them with the distance travelled to have the neutron's motion.
Reflection 4: Simulation phyiscs and investigation (B)¶
(Highlight something you feel you did particularily well in terms of the context of your simulation, the physical modelling that you did or the investigation you performed. Discuss an aspect of these dimensions of your project that would benefit the most from further effort.)
For this project, the simulation that I think I did well is how I generated uniform random motion of each neutron using the method I mentioned in Reflection 3. For a more precise simulation, I think I will need bigger trials since 250 trials might be not enough and cause the mean k factor to deviate from the true/expected value.
Reflection 5: Effectiveness of your communication¶
(Highlight something you feel you did particularily well in your visualizations or written communication. Discuss an aspect of your visualizations or written communication that would benefit the most from further effort.)
For this project, I learned from my previous projects to add very specific comments and notes in each important step of code so that it is less confusing for me and TA to revise it. On the other hand, my visualizations can be improved by adding more details and labels such as the corresponding a (side-length) of each volume value in the body plot.