<jupyter_start><jupyter_text>Lab 5 - Modified 1-Group CriticalityIn this lab we will do a modified 1-group criticality calculation for a homogeneous thermal reactor comprising a sphere of U-235 with light water.  The composition is specified (mass density of U-235) and the goal will be to find the critical radius of the sphere. The calculation will be repeated to take into account the effect of a light water reflector.In order to compare the diffusion theory result with that obtained using a neutron transport simulation package, we repeat the calculation using OpenMC.  We will also find the critical radius in the case where the reflector material is changed from light water to graphite.In all cases it will be assumed that the thermal medium is maintained at 20${\circ}C.<jupyter_code>import numpy as np

# 1-Group solution

# Parameters
N_M = 0.03343; # atoms/b-cm, atom density of water moderator at 20 C

rho_F = 0.0145; # g/cm**3, mass density of the fuel (given)
N_A = 0.6022; # atoms-cm**2/mol-b, Avagadro's number in weird units
M_F = 235.0; # g/mol, gram atomic weight of fuel (approximate)
N_F = rho_F*N_A/M_F; # atoms/b-cm, atom density of the fuel

sig_a_F = 687; # b, 0.0253-eV micro absorption cros section for U-235
g_a_F = 0.9780; # non 1-v factor for U-235 absorption at 20 C
sig_a_M = 0.664; # b, 0.0253-ev micro absorption cross section for light water
eta_F = 2.065; # reproduction factor for U-235 at 20 C
tau_T = 27; # cm**2, fermi age for water at 20 C
L_T_M_sq = 8.1; # cm**2, thermal diffusion area for water at 20 C

Z = N_F*sig_a_F*g_a_F/(N_M*sig_a_M); # modified 1-group (convenience) parameter. 
#note constant corrections to average thermal absorption cross sections cancel out.
f = Z/(Z+1.); # thermal utilization
L_T_sq = L_T_M_sq/(Z+1.); # thermal diffusion area accounting for presence of fuel
D_M = 0.16; # cm, diffusion area of moderator at 20 C
d = 2.13*D_M; # cm, extrapolation distance

k_inf = eta_F*f; # inifinite multiplication factor
B_sq = (k_inf - 1.)/(L_T_sq + tau_T); # 1/cm**2, critical buckling
R_extrap = np.pi/np.sqrt(B_sq);


R = R_extrap - d; # cm, critical radius

print(f'Critical Radius from modified 1-group theory: %5.4f cm'%R);


# Reflector Savings
delta = 7.2 + 0.10*((L_T_sq+tau_T)-40.0); # correlation for light water moderated/reflected cores

print(f'Estimated reflector savings: %5.4f cm'%delta);

R_ref = R - delta;
print(f'Estimated reflected critical radius: %5.4f cm'%R_ref);<jupyter_output>Estimated reflector savings: 6.2812 cm
Estimated reflected critical radius: 50.5330 cm<jupyter_text>OpenMC Model<jupyter_code>import openmc
# create a model for the critical radius search
rho_M = 1.0
rho_tot = rho_M+rho_F;

wo_F = rho_F/rho_tot;
wo_H = (2./18.)*rho_M/rho_tot;
wo_O = (16./18.)*rho_M/rho_tot;


# create the model. "R" will be the parametric variable
def build_model(R):
    # U-235 and light water sphere with given fuel composition
    fuel = openmc.Material(name='fuel');
    fuel.add_nuclide('U235',wo_F,'wo');
    fuel.add_nuclide('H1',wo_H,'wo'); # 2 H1's per H2O molecule
    fuel.add_nuclide('O16',wo_O,'wo');# 1 O16's per H2O molecule
    fuel.add_s_alpha_beta('c_H_in_H2O'); # S(alpha,beta) for water.
    fuel.set_density('g/cm3',1.0+rho_F); # approximate mass density of mixture
    
    materials = openmc.Materials([fuel]);
    
    fuel_sphere = openmc.Sphere(r=R,boundary_type='vacuum');
    
    fuel_cell = openmc.Cell();
    fuel_cell.fill = fuel;
    fuel_cell.region = -fuel_sphere;
    
    root_universe = openmc.Universe();
    root_universe.add_cells([fuel_cell]);
    
    geometry = openmc.Geometry(root_universe);
    
    settings = openmc.Settings()
    settings.batches = 200;
    settings.inactive = 50;
    settings.particles = 10000;
    
    bounds = [-R,-R,-R,R,R,R];
    uniform_dist = openmc.stats.Box(bounds[:3],bounds[3:],
                                   only_fissionable=True);
    settings.source = openmc.source.Source(space=uniform_dist);
    
    settings.output = {'tallies':False};
    
    model = openmc.model.Model(geometry,materials,settings);
    
    return model;
crit_R, guesses, keffs = openmc.search_for_keff(build_model,
                                                bracket=[40.,65.],
                                                tol=1e-3,
                                                print_iterations=True);
print(f'Critical Radius: %5.4f cm' % crit_R);<jupyter_output>Iteration: 1; Guess of 4.00e+01 produced a keff of 0.93677 +/- 0.00065
Iteration: 2; Guess of 6.50e+01 produced a keff of 1.02429 +/- 0.00056
Iteration: 3; Guess of 5.25e+01 produced a keff of 0.99404 +/- 0.00065
Iteration: 4; Guess of 5.88e+01 produced a keff of 1.01344 +/- 0.00074
Iteration: 5; Guess of 5.56e+01 produced a keff of 1.00409 +/- 0.00068
Iteration: 6; Guess of 5.41e+01 produced a keff of 0.99894 +/- 0.00071
Iteration: 7; Guess of 5.48e+01 produced a keff of 1.00148 +/- 0.00063
Iteration: 8; Guess of 5.45e+01 produced a keff of 1.00045 +/- 0.00067
Iteration: 9; Guess of 5.43e+01 produced a keff of 0.99993 +/- 0.00070
Iteration: 10; Guess of 5.44e+01 produced a keff of 1.00016 +/- 0.00071
Iteration: 11; Guess of 5.43e+01 produced a keff of 1.00015 +/- 0.00067
Critical Radius: 54.3066 cm<jupyter_text>This isn't too far from the MCNP result.  If your hairs stand on end when you look at the way that I handled the fuel composition then you'll understand why I'm not super confident about this result, but it doesn't look bad to me.Update:  Actually, I probably did a better job of defining the material in this than what I did in the MCNP input (if you can believe that).  The first time around (i.e. when I wrote the note above) I had neglected to include the S(alpha,beta) thermal scattering treatment in the fuel material.  That resulted in ~2 cm difference in the computed critical radius.  I declare that to be non-trivial.Lesson learned: pay attention to the details!!<jupyter_code># model with light water reflector
rho_M = 1.0
rho_tot = rho_M+rho_F;

wo_F = rho_F/rho_tot;
wo_H = (2./18.)*rho_M/rho_tot;
wo_O = (16./18.)*rho_M/rho_tot;


# create the model. "R" will be the parametric variable
def water_reflected_model(R):
    # U-235 and light water sphere with given fuel composition
    fuel = openmc.Material(name='fuel');
    fuel.add_nuclide('U235',wo_F,'wo');
    fuel.add_nuclide('H1',wo_H,'wo'); # 2 H1's per H2O molecule
    fuel.add_nuclide('O16',wo_O,'wo');# 1 O16's per H2O molecule
    fuel.add_s_alpha_beta('c_H_in_H2O'); # S(alpha,beta) for water.
    fuel.set_density('g/cm3',1.0+rho_F); # approximate mass density of mixture
    
    lw_ref = openmc.Material(name='light water reflector');
    lw_ref.add_nuclide('H1',2.,'ao');
    lw_ref.add_nuclide('O16',1.,'ao');
    lw_ref.add_s_alpha_beta('c_H_in_H2O');
    lw_ref.set_density('g/cm3',1.0);
    
    
    materials = openmc.Materials([fuel,lw_ref]);
    
    fuel_sphere = openmc.Sphere(r=R);
    outer_sphere = openmc.Sphere(r=5.*R,boundary_type='vacuum');
    # assume the reflector is 5x the radius of the sphere.  That should be enough 
    # for an (effectively) infinite reflector
    
    fuel_cell = openmc.Cell();
    fuel_cell.fill = fuel;
    fuel_cell.region = -fuel_sphere;
    
    ref_cell = openmc.Cell();
    ref_cell.fill = lw_ref;
    ref_cell.region = +fuel_sphere & -outer_sphere;
    
    root_universe = openmc.Universe();
    root_universe.add_cells([fuel_cell,ref_cell]);
    
    geometry = openmc.Geometry(root_universe);
    
    settings = openmc.Settings()
    settings.batches = 200;
    settings.inactive = 50;
    settings.particles = 10000;
    
    bounds = [-R,-R,-R,R,R,R];
    uniform_dist = openmc.stats.Box(bounds[:3],bounds[3:],
                                   only_fissionable=True);
    settings.source = openmc.source.Source(space=uniform_dist);
    
    settings.output = {'tallies':False};
    
    model = openmc.model.Model(geometry,materials,settings);
    
    return model;
crit_R_reflected, guesses, keffs = openmc.search_for_keff(water_reflected_model,
                                                bracket=[40.,65.],
                                                tol=1e-3,
                                                print_iterations=True);
print(f'Critical Radius: %5.4f cm' % crit_R_reflected);<jupyter_output>Iteration: 1; Guess of 4.00e+01 produced a keff of 0.95613 +/- 0.00070
Iteration: 2; Guess of 6.50e+01 produced a keff of 1.03124 +/- 0.00062
Iteration: 3; Guess of 5.25e+01 produced a keff of 1.00521 +/- 0.00056
Iteration: 4; Guess of 4.62e+01 produced a keff of 0.98413 +/- 0.00067
Iteration: 5; Guess of 4.94e+01 produced a keff of 0.99470 +/- 0.00066
Iteration: 6; Guess of 5.09e+01 produced a keff of 1.00022 +/- 0.00060
Iteration: 7; Guess of 5.02e+01 produced a keff of 0.99830 +/- 0.00059
Iteration: 8; Guess of 5.05e+01 produced a keff of 0.99940 +/- 0.00067
Iteration: 9; Guess of 5.07e+01 produced a keff of 0.99915 +/- 0.00064
Iteration: 10; Guess of 5.08e+01 produced a keff of 0.99951 +/- 0.00069
Iteration: 11; Guess of 5.09e+01 produced a keff of 0.99846 +/- 0.00069
Critical Radius: 50.8887 cm<jupyter_text>This is, again, similar to my previously computed MCNP result.  Note that the reflector savings is not as big as predicted by the correlation we used for the "analytic" method.<jupyter_code># model with graphite reflector
rho_M = 1.0
rho_tot = rho_M+rho_F;

wo_F = rho_F/rho_tot;
wo_H = (2./18.)*rho_M/rho_tot;
wo_O = (16./18.)*rho_M/rho_tot;


# create the model. "R" will be the parametric variable
def graphite_reflected_model(R):
    # U-235 and light water sphere with given fuel composition
    fuel = openmc.Material(name='fuel');
    fuel.add_nuclide('U235',wo_F,'wo');
    fuel.add_nuclide('H1',wo_H,'wo'); # 2 H1's per H2O molecule
    fuel.add_nuclide('O16',wo_O,'wo');# 1 O16's per H2O molecule
    fuel.add_s_alpha_beta('c_H_in_H2O'); # S(alpha,beta) for water.
    fuel.set_density('g/cm3',1.0+rho_F); # approximate mass density of mixture
    
    grph_ref = openmc.Material(name='graphite reflector');
    grph_ref.add_nuclide('C0',1.,'ao');
    grph_ref.add_s_alpha_beta('c_Graphite');
    grph_ref.set_density('g/cm3',1.7);
    
    
    materials = openmc.Materials([fuel,grph_ref]);
    
    fuel_sphere = openmc.Sphere(r=R);
    outer_sphere = openmc.Sphere(r=5.*R,boundary_type='vacuum');
    # assume the reflector is 5x the radius of the sphere.  That should be enough 
    # for an (effectively) infinite reflector
    
    fuel_cell = openmc.Cell();
    fuel_cell.fill = fuel;
    fuel_cell.region = -fuel_sphere;
    
    ref_cell = openmc.Cell();
    ref_cell.fill = grph_ref;
    ref_cell.region = +fuel_sphere & -outer_sphere;
    
    root_universe = openmc.Universe();
    root_universe.add_cells([fuel_cell,ref_cell]);
    
    geometry = openmc.Geometry(root_universe);
    
    settings = openmc.Settings()
    settings.batches = 200;
    settings.inactive = 50;
    settings.particles = 10000;
    
    bounds = [-R,-R,-R,R,R,R];
    uniform_dist = openmc.stats.Box(bounds[:3],bounds[3:],
                                   only_fissionable=True);
    settings.source = openmc.source.Source(space=uniform_dist);
    
    settings.output = {'tallies':False};
    
    model = openmc.model.Model(geometry,materials,settings);
    
    return model;
crit_R_grph_ref, guesses, keffs = openmc.search_for_keff(graphite_reflected_model,
                                                bracket=[40.,65.],
                                                tol=1e-3,
                                                print_iterations=True);
print(f'Critical Radius: %5.4f cm' % crit_R_grph_ref);<jupyter_output>Iteration: 1; Guess of 4.00e+01 produced a keff of 0.99512 +/- 0.00063
Iteration: 2; Guess of 6.50e+01 produced a keff of 1.04948 +/- 0.00059
Iteration: 3; Guess of 5.25e+01 produced a keff of 1.02955 +/- 0.00061
Iteration: 4; Guess of 4.62e+01 produced a keff of 1.01607 +/- 0.00058
Iteration: 5; Guess of 4.31e+01 produced a keff of 1.00605 +/- 0.00068
Iteration: 6; Guess of 4.16e+01 produced a keff of 1.00154 +/- 0.00066
Iteration: 7; Guess of 4.08e+01 produced a keff of 0.99844 +/- 0.00063
Iteration: 8; Guess of 4.12e+01 produced a keff of 1.00079 +/- 0.00067
Iteration: 9; Guess of 4.10e+01 produced a keff of 1.00032 +/- 0.00060
Iteration: 10; Guess of 4.09e+01 produced a keff of 0.99926 +/- 0.00064
Iteration: 11; Guess of 4.09e+01 produced a keff of 0.99856 +/- 0.00063
Iteration: 12; Guess of 4.10e+01 produced a keff of 0.99947 +/- 0.00062
Critical Radius: 40.9521 cm
{'fraction_non_alphanumeric': 0.11548577349167131, 'fraction_numerical': 0.07404160357057464, 'mean_word_length': 3.1426213271706835, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '79774', 'license': 'MIT', 'n_tokens_mistral': 5170, 'n_tokens_neox': 4397, 'n_words': 1271, 'path': 'lab5/lab5.ipynb', 'repo_name': 'stu314159/er362_omc', 'size': '12547'}