<jupyter_start><jupyter_text>Plutonium MakerThe purpose of this notebook is to investigate the performance of various pin-cell configurations in producing weapons grade plutonium.The sorts of questions that I want to be able to answer:1. What are the trade-offs in moderator selection2. What should be the fuel pitch, diameter and/or pitch to diameter ratio?3. How long should the fuel be left in the reactor for maximum production of plutonium? 4.  I expect the answer to 3 to depend on the power density and the neutron energy spectrum.  I assume you should just operate at the maximum achievable power density (with respect to heat transfer requirements) but how does the neutron energy spectrum change things?5. If I wanted to make Pu-238, how would the reactor core design change. (fuel, moderator, coolant, geometric arrangement, etc...)After some initial investigations, I can investigate the available design information for early plutonium production reactors.  As much as is possible, I could try to observe how plutonium production reactor core design has changed over the years and try to infer the engineering rationale for the evolution. In short: I hope to build a *quantitative* intuition regarding these considerations.<jupyter_code>import openmc
import numpy as np

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Materials<jupyter_code>fuel = openmc.Material(name="uo2")
fuel.add_element("U", 1, percent_type="ao", enrichment=0.7)
#fuel.add_element("O", 2)
fuel.set_density("g/cc", 18.4)

clad = openmc.Material(name='clad');
clad.add_element("Zr",1);
clad.set_density('g/cc',6.0);

alum = openmc.Material(name='alum');
alum.add_element('Al',1);
alum.set_density('g/cc',2.7); #according to Google...

water = openmc.Material(name='water');
water.add_element('O',1);
water.add_element('H',2)
water.set_density('g/cc',0.712); # higher temperature water...
water.add_s_alpha_beta('c_H_in_H2O');

graphite = openmc.Material(name='graphite');
graphite.add_element('C',1);
graphite.set_density('g/cc',1.7);
graphite.add_s_alpha_beta('c_Graphite');

heavy_w = openmc.Material(name='heavy_w');
heavy_w.add_element('O',1);
heavy_w.add_nuclide('H2',2);
heavy_w.set_density('g/cc',1.1);
heavy_w.add_s_alpha_beta('c_D_in_D2O');

materials = openmc.Materials([fuel,clad,water,alum,graphite,heavy_w]);
materials.export_to_xml();<jupyter_output><empty_output><jupyter_text>GeometryI will make a pin cell.  In order to (roughly) model a hanford pin cell I will use a uranium metal slug coated with aluminum in a graphite reflector cooled with water.<jupyter_code>h_core = 300.; # all dimensions in cm 
h_fuel = 200.;

r_fuel = 2.0;
r_pin = 2.1;
r_cool = 2.25;

P_D = 6; 

pitch = P_D*r_cool*2.; # set P/D relative to the coolant channel diameter
# P_D=2 is too tight; k_eff = 0.93
# on later iterations add temperature effects to the fuel to see how important they are.

fuel_cyl = openmc.model.RightCircularCylinder([0.,0.,-h_fuel/2.],
                                             h_fuel, r_fuel);
fuel.volume = np.pi*(r_fuel**2)*h_fuel;

pin_cyl = openmc.model.RightCircularCylinder([0.,0.,-(h_fuel+(r_pin-r_fuel))/2.],
                                            h_fuel+(r_pin-r_fuel)*2.,r_pin);

cool_ch_cyl = openmc.model.RightCircularCylinder([0.,0.,-h_core/2.+5.],
                                                h_core-10.,r_cool);
# note: I added some margin between the coolant cell upper and lower boundaries and the core.
# having those surfaces coincident resulted in lost particles.

core_cell = openmc.model.RectangularParallelepiped(-pitch/2.,pitch/2.,
                                                  -pitch/2.,pitch/2.,
                                                  -h_core/2.,h_core/2.,
                                                  boundary_type='reflective');

fuel_cell = openmc.Cell();
fuel_cell.region = -fuel_cyl;
fuel_cell.fill = fuel;
fuel_cell.temperature = 800; #K, reminder to fill this in later

clad_cell = openmc.Cell();
clad_cell.region = +fuel_cyl & -pin_cyl
clad_cell.fill = clad;

cool_cell = openmc.Cell();
cool_cell.region = +pin_cyl & -cool_ch_cyl;
cool_cell.fill = water;

mod_cell = openmc.Cell();
mod_cell.region = +cool_ch_cyl & -core_cell;
mod_cell.fill = heavy_w;
mod_cell.temperature = 300; # K, just a guess



root_univ = openmc.Universe();
root_univ.add_cells([fuel_cell,clad_cell,cool_cell,mod_cell]);

geometry = openmc.Geometry();
geometry.root_universe = root_univ;

geometry.export_to_xml();

color = {}
color[fuel]='yellow';
color[alum]='gray';
color[water]='blue';
color[graphite]='black';
color[heavy_w] = 'cyan';

#root_univ.plot(width=(pitch,pitch),color_by='material',colors=color);
root_univ.plot(width=(3*r_cool,3*r_cool),color_by='material',colors=color);<jupyter_output><empty_output><jupyter_text>SettingsInitially we will just find what the eigenvalue is.  Since this will model a reactor that is infinite, the materials will need to be tweaked so that the eigenvalue is greater than 1.<jupyter_code>settings = openmc.Settings();
settings.run_mode = 'eigenvalue';
settings.particles = 10000;
settings.batches = 150;
settings.inactive = 50;

box = openmc.stats.Box(lower_left = (-r_fuel,-r_fuel,-h_fuel/2.),
                      upper_right = (r_fuel,r_fuel,h_fuel/2.),
                      only_fissionable=True);
src = openmc.Source(space=box);

settings.source = src;

settings.temperature['method']='interpolation';

settings.export_to_xml();

openmc.run();<jupyter_output>%%%%%%%%%%%%%%%
                           %%%%%%%%%%%%%%%%%%%%%%%%
                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
                                    %%%%%%%%%%%%%%%%%%%%%%%%
                                     %%%%%%%%%%%%%%%%%%%%%%%%
                 ###############      %%%%%%%%%%%%%%%%%%%%%%%%
                ##################     %%%%%%%%%%%%%%%%%%%%%%%
                ###################     %%%%%%%%%%%%%%%%%%%%%%%
                ####################     %%%%%%%%%%%%%%%%%%%%%%
                #####################     %%%%%%%%%%%%%%%%%%%%%
                ######################     %%%%%%%%%%%%%%%%%%%%
                #######################     %%%%%%%%%%%%%%%%%%
                 #######################     %%%%%%%%%%%%%%%%%
                 #####################[...]<jupyter_text>Depletion CalculationOkay - this is a start.  I fully expect that the reactor, as modeled, will not remain critical once equilibrium Xe-135 and Sm-149 are established in the core.  The negative reactivity from the doppler broadening of the fuel didn't turn out to be such a big deal.<jupyter_code>import openmc.deplete

operator = openmc.deplete.Operator(geometry,settings,"chain_casl_pwr.xml");
power = 1e4; # W, this is a guess

operator.heavy_metal

max_step = 2*operator.heavy_metal/power*1e3
print(f'Maximum recommended depletion step %5.3f days' % max_step);

days = 24*3600; # seconds in a day
time_steps = [0.1*days,0.1*days,0.2*days,0.3*days,0.5*days,
              1.0*days,10.0*days, 100.0*days, 100.0*days, 200.*days,
             200.*days, 200.*days, 500.*days, 500.*days];

integrator = openmc.deplete.PredictorIntegrator(operator,time_steps,power=power);
#integrator = openmc.deplete.CELIIntegrator(operator,time_steps,power=power);
#integrator = openmc.deplete.SICELIIntegrator(operator,time_steps,power=power);
integrator.integrate()
results = openmc.deplete.ResultsList.from_hdf5("./depletion_results.h5")
time, k = results.get_eigenvalue()
time /= (24*3600);
from matplotlib import pyplot
pyplot.errorbar(time,k[:,0],yerr=k[:,1]);
pyplot.xlabel('Time [d]');
pyplot.ylabel('$k_{eff}\pm \sigma$');
pyplot.title("$k_{eff}$ vs time");
pyplot.grid();<jupyter_output><empty_output><jupyter_text>obviously, these results suck.  A better analysis would use more particles at least and probably smaller time steps.<jupyter_code>_time, pu239 = results.get_atoms("1","Pu239");
_time, pu240 = results.get_atoms("1","Pu240");
_time, pu241 = results.get_atoms("1","Pu241");
pyplot.semilogy(time,pu239,label="Pu239");
pyplot.semilogy(time,pu240,label="Pu240");
pyplot.semilogy(time,pu241,label="Pu241");
pyplot.grid();
pyplot.xlabel('Time [d]');
pyplot.ylabel('Number of atoms');
pyplot.legend();

mass_pu239 = pu239*(239./6.022e23);
mass_pu240 = pu240*(240./6.022e23);
mass_pu241 = pu241*(241./6.022e23);

mass_pu = mass_pu239 + mass_pu240 + mass_pu241
pu_239_ratio = mass_pu239/mass_pu
print(pu_239_ratio)
_time, xe135 = results.get_atoms("1","Xe135")
_time, sm149 = results.get_atoms("1","Sm149")

pyplot.semilogy(time,xe135,label="Xe135");
pyplot.semilogy(time,sm149,label="Sm149");
pyplot.grid();
pyplot.xlabel('Time [d]');
pyplot.ylabel('Number of atoms');
pyplot.legend();

results_list = openmc.deplete.ResultsList.<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18236877523553163, 'fraction_numerical': 0.03297442799461642, 'mean_word_length': 4.006737787759685, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 24, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 4}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '83326', 'license': 'MIT', 'n_tokens_mistral': 3133, 'n_tokens_neox': 2772, 'n_words': 643, 'path': 'lab8/pu_maker/hw_version/pu_breeder.ipynb', 'repo_name': 'stu314159/er362_omc', 'size': '8916'}