I have a model of a CANDU reactor that I am trying to deplete and then perform a boron criticality search before running another depletion simulation, and when refueling begins it will run depletion->refueling->criticality search using adjuster rods. I am having issues with the script to perform this. I have a successful boron criticality search script, along with the depletion, so I thought it would be as easy as putting them in a loop, but I am having issues with properly updating the boron concentration after a depletion is ran. Below is the entire script of a simplified model that I am using as a replacement to accelerate the trial and error process.
time_steps = [5/24] * 8
current_time = 0.0
for i, dt in enumerate(time_steps, start=1):
current_time += dt
step_id = f"{current_time:.10f}_days" #Assigns the step id to the number of days
model = openmc.Model(
geometry=geometry,
materials=materials,
settings=settings,
# tallies=tallies_reactionR_file
)
operator = openmc.deplete.CoupledOperator(model, "chain.xml", normalization_mode="energy-deposition")
integrator = openmc.deplete.PredictorIntegrator(operator, [time_steps[i-1]], power, timestep_units='d', solver='cram16').integrate()
rename_and_move_openmc_outputs(step_id)
filename = f"statepoint.20_{step_id}.h5"
# construct the statepoint path to call in later statepoints
statepoint_path = os.path.join(parent_folder_name, step_id, filename)
results = openmc.deplete.Results(f'{parent_folder_name}/{step_id}/depletion_results_{step_id}.h5')
materials = results.export_to_materials(1)
#Arguments for the boron criticality function
Boron_function_args = {
#Dimensions in cm
'r_fuel': r_fuel,
'r_clad': r_clad,
'pitch': pitch,
'pressure_tube_ir': pressure_tube_ir,
'pressure_tube_or': pressure_tube_or,
'calandria_ir': calandria_ir,
'calandria_or': calandria_or,
'fuel_bundle_length': fuel_bundle_length,
'sheath_length': sheath_length,
'r_coolant': r_coolant,
'Calandria_vessel_ir': Calandria_vessel_ir,
'Calandria_vessel_or': Calandria_vessel_or,
'Calandria_vessel_length': Calandria_vessel_length,
'reflective_distance_x': reflective_distance_x,
'reflective_distance_y': reflective_distance_y,
'void_z': void_z,
#Temperatures in Kelvin
'moderator_temp': moderator_temp,
'fuel_temp': fuel_temp,
'coolant_temp': coolant_temp,
'clad_temp': clad_temp,
'pressure_tube_temp': pressure_tube_temp,
'anulus_temp': anulus_temp,
'calandria_tube_temp': calandria_tube_temp,
'calandria_vessel_temp': calandria_vessel_temp,
'adjuster_temp': adjuster_temp,
'LCZ_temp': LCZ_temp,
#Densities in g/cm3
'moderator_dens': moderator_dens,
'fuel_dens': fuel_dens,
'coolant_dens': coolant_dens,
'clad_dens': clad_dens,
'pressure_tube_dens': pressure_tube_dens,
'anulus_dens': anulus_dens,
'calandria_tube_dens': calandria_tube_dens,
'calandria_vessel_dens': calandria_vessel_dens,
'adjuster_dens': adjuster_dens,
'LCZ_dens': LCZ_dens,
#D2O Atomic Concentraion in homogenized adjuster rods between 0 and 1
'D2OC': D2OC,
#Settings
'batches': batches,
#Need at least 100 entries removed, so 10 batches*10 gens, or 100 batches*1 gen
'inactive_batches': inactive_batches,
'generations': generations,
'particles_boron': particles_boron,
#Misc
'boron_enrichment': boron_enrichment,
'fuel_materials' : [],
'fuel_array' : fuel_array
}
#Parametized model for boron criticality search
def Boron_build_model(ppm_Boron, **kwargs):
openmc.reset_auto_ids() # clears ID counters and registries
r_fuel = kwargs["r_fuel"]
r_clad = kwargs["r_clad"]
pitch = kwargs["pitch"]
pressure_tube_ir = kwargs["pressure_tube_ir"]
pressure_tube_or = kwargs["pressure_tube_or"]
calandria_ir = kwargs["calandria_ir"]
calandria_or = kwargs["calandria_or"]
fuel_bundle_length = kwargs["fuel_bundle_length"]
sheath_length = kwargs["sheath_length"]
r_coolant = kwargs["r_coolant"]
Calandria_vessel_ir = kwargs["Calandria_vessel_ir"]
Calandria_vessel_or = kwargs["Calandria_vessel_or"]
Calandria_vessel_length = kwargs["Calandria_vessel_length"]
reflective_distance_x = kwargs["reflective_distance_x"]
reflective_distance_y = kwargs["reflective_distance_y"]
void_z = kwargs["void_z"]
#Temperatures in Kelvin
moderator_temp = kwargs["moderator_temp"]
fuel_temp = kwargs["fuel_temp"]
coolant_temp = kwargs["coolant_temp"]
clad_temp = kwargs["clad_temp"]
pressure_tube_temp = kwargs["pressure_tube_temp"]
anulus_temp = kwargs["anulus_temp"]
calandria_tube_temp = kwargs["calandria_tube_temp"]
calandria_vessel_temp = kwargs["calandria_vessel_temp"]
adjuster_temp = kwargs["adjuster_temp"]
LCZ_temp = kwargs["LCZ_temp"]
#Densities in g/cm3
moderator_dens = kwargs["moderator_dens"]
fuel_dens = kwargs["fuel_dens"]
coolant_dens = kwargs["coolant_dens"]
clad_dens = kwargs["clad_dens"]
pressure_tube_dens = kwargs["pressure_tube_dens"]
anulus_dens = kwargs["anulus_dens"]
calandria_tube_dens = kwargs["calandria_tube_dens"]
calandria_vessel_dens = kwargs["calandria_vessel_dens"]
adjuster_dens = kwargs["adjuster_dens"]
LCZ_dens = kwargs["LCZ_dens"]
#D2O Atomic Concentraion in homogenized adjuster rods between 0 and 1
D2OC = kwargs["D2OC"]
#Settings
batches = kwargs["batches"]
#Need at least 100 entries removed, so 10 batches*10 gens, or 100 batches*1 gen
inactive_batches = kwargs["inactive_batches"]
generations = kwargs["generations"]
particles_boron = kwargs["particles_boron"]
#Misc
boron_enrichment = kwargs["boron_enrichment"]
#arrays used for fuel
fuel_materials = []
fuel_array = kwargs["fuel_array"]
#Materials Definition and Export
Sheath = openmc.Material(1,"Sheath")
Sheath.add_nuclide('Zr90', 0.981)
Sheath.add_element('Sn', 0.014)
Sheath.add_element('Fe', 0.002)
Sheath.add_element('Cr', 0.001)
Sheath.add_element('O', 0.002)
Sheath.set_density('g/cm3', clad_dens)
Calandria_tube = openmc.Material(2,'Calandria Tube')
Calandria_tube.add_nuclide('Zr90', 0.9823)
Calandria_tube.add_element('Sn', 0.014)
Calandria_tube.add_element('Fe', 0.001)
Calandria_tube.add_element('Cr', 0.001)
Calandria_tube.add_element('O', 0.0012)
Calandria_tube.add_element('Ni', 0.0005)
Calandria_tube.set_density('g/cm3', calandria_tube_dens)
Pressure_tube = openmc.Material(3,"Pressure Tube")
Pressure_tube.add_nuclide('Zr90', 0.975)
Pressure_tube.add_element('Nb', 0.025)
Pressure_tube.set_density('g/cm3', pressure_tube_dens)
Moderator = openmc.Material(4,"Moderator")
Moderator.add_nuclide('H2', 1.995)
Moderator.add_nuclide('H1', 0.005)
Moderator.add_nuclide('O16', 1.0)
Moderator.set_density('g/cm3', moderator_dens)
Moderator.add_element('B', ppm_Boron, enrichment=boron_enrichment, enrichment_target='B10')
Moderator.add_s_alpha_beta('c_D_in_D2O')
Coolant = openmc.Material(5, "Coolant")
Coolant.add_nuclide('H2', 1.98)
Coolant.add_nuclide('H1', 0.02)
Coolant.add_nuclide('O16', 1.0)
Coolant.set_density('g/cm3', coolant_dens)
Coolant.add_s_alpha_beta('c_D_in_D2O')
CO2 = openmc.Material(6, "CO2")
CO2.add_element('C', 1.0)
CO2.add_element('O', 2.0)
CO2.set_density('g/cm3', anulus_dens)
Calandria_vessel = openmc.Material(7, "Calandria vessel")
Calandria_vessel.add_element('C', 0.0003)
Calandria_vessel.add_element('Mn', 0.02)
Calandria_vessel.add_element('P', 0.00045)
Calandria_vessel.add_element('S', 0.0003)
Calandria_vessel.add_element('Si', 0.0075)
Calandria_vessel.add_element('Cr', 0.19)
Calandria_vessel.add_element('Ni', 0.10)
Calandria_vessel.add_element('N', 0.0010)
Calandria_vessel.add_element('Fe', 0.68045)
Calandria_vessel.set_density('g/cm3', calandria_vessel_dens)
Air = openmc.Material(8, "Air")
Air.add_element('N', 0.7812)
Air.add_element('O', 0.2095)
Air.add_element('Ar', 0.0093)
Air.set_density('g/cm3', 0.001225)
AdjusterA = openmc.Material(9, "AdjusterA")
AdjusterA.add_element('C', 0.0003)
AdjusterA.add_element('Mn', 0.02)
AdjusterA.add_element('P', 0.00045)
AdjusterA.add_element('S', 0.0003)
AdjusterA.add_element('Si', 0.0075)
AdjusterA.add_element('Cr', 0.19)
AdjusterA.add_element('Ni', 0.10)
AdjusterA.add_element('N', 0.0010)
AdjusterA.add_element('Fe', 0.68045)
AdjusterA.set_density('g/cm3', adjuster_dens)
#Mixing Homogenized Adjuster composition
ModeratorA = openmc.Material(10,"ModeratorA")
ModeratorA.add_nuclide('H2', 1.995)
ModeratorA.add_nuclide('H1', 0.005)
ModeratorA.add_nuclide('O16', 1.0)
ModeratorA.set_density('g/cm3', moderator_dens)
Adjuster = openmc.Material.mix_materials([AdjusterA, ModeratorA], [1-D2OC, D2OC], 'ao')
Adjuster.add_s_alpha_beta('c_D_in_D2O', fraction=D2OC)
Helium = openmc.Material(12, "Helium")
Helium.add_element('Helium', 2.0)
Helium.set_density('g/cm3', 0.000037239)
Light_water = openmc.Material(13, "Light Water")
Light_water.add_nuclide('H1', 2.0)
Light_water.add_nuclide('O16', 1.0)
Light_water.set_density('g/cm3', LCZ_dens)
Light_water.add_s_alpha_beta('c_H_in_H2O')
for i, code in enumerate(fuel_array, start=14):
# Create the variable name dynamically by appending the code
exec(f"Fuel_{code} = openmc.Material({i}, 'Fuel_{code}')")
#Execute the OpenMC commands for the dynamically created variable
exec(f"Fuel_{code}.add_element('U', 1.0, enrichment=0.71)")
exec(f"Fuel_{code}.add_nuclide('O16', 2.0)")
exec(f"Fuel_{code}.set_density('g/cm3', fuel_dens)")
exec(f"Fuel_{code}.volume = fuel_volume")
exec(f"fuel_materials.append(Fuel_{code})")
materials = openmc.Materials([Sheath, Calandria_tube, Pressure_tube, Moderator, Coolant, CO2, Calandria_vessel, Air, Adjuster, Helium, Light_water, AdjusterA] + fuel_materials)
materials.cross_sections = '/home/smrcenter/openmc/endfb80/cross_sections.xml'
materials.export_to_xml()
#Making planes and cylinders for creating fuel geometry
surf_fuel = openmc.ZCylinder(r=r_fuel)
back_surf_fuel = openmc.ZPlane(-fuel_bundle_length/2)
front_surf_fuel = openmc.ZPlane(fuel_bundle_length/2)
back_surf_clad = openmc.ZPlane(-sheath_length/2)
front_surf_clad = openmc.ZPlane(sheath_length/2)
clad_outside_surface = openmc.ZCylinder(r=r_clad)
coolant_inner_surf = openmc.ZCylinder(r=r_coolant)
pt_inner = openmc.ZCylinder(r=pressure_tube_ir)
pt_outer = openmc.ZCylinder(r=pressure_tube_or)
calandria_inner = openmc.ZCylinder(r=calandria_ir)
calandria_outer = openmc.ZCylinder(r=calandria_or)
#Pitch box
left = openmc.XPlane(-pitch/2)
right = openmc.XPlane(pitch/2)
bottom = openmc.YPlane(-pitch/2)
top = openmc.YPlane(pitch/2)
back = openmc.ZPlane(-sheath_length/2)
front = openmc.ZPlane(sheath_length/2)
pitch_box = +left & -right & +bottom & -top & +back & -front
#Making Fuel cells
for i, code in enumerate(fuel_array):
exec(f"fuel_cell_{code} = openmc.Cell(fill=Fuel_{code}, region=-surf_fuel & +back_surf_fuel & -front_surf_fuel)")
exec(f"fuel_cell_{code}.temperature = fuel_temp")
# Create clad side cell for each material
exec(f"clad_side_cell_{code} = openmc.Cell(fill=Sheath, region= +surf_fuel & -clad_outside_surface & +back_surf_fuel & -front_surf_fuel)")
exec(f"clad_side_cell_{code}.temperature = clad_temp")
# Create clad back cell for each material
exec(f"clad_back_cell_{code} = openmc.Cell(fill=Sheath, region=-clad_outside_surface & -back_surf_fuel & +back_surf_clad)")
exec(f"clad_back_cell_{code}.temperature = clad_temp")
# Create clad front cell for each material
exec(f"clad_front_cell_{code} = openmc.Cell(fill=Sheath, region=-clad_outside_surface & +front_surf_fuel & -front_surf_clad)")
exec(f"clad_front_cell_{code}.temperature = clad_temp")
exec(f"coolant_cell_inner_{code} = openmc.Cell(fill=Coolant, region=-coolant_inner_surf & +back_surf_clad & -front_surf_clad)")
exec(f"coolant_cell_inner_{code}.temperature = coolant_temp")
exec(f"coolant_cell_outer_{code} = openmc.Cell(fill=Coolant, region=+coolant_inner_surf & -pt_inner & +back_surf_clad & -front_surf_clad)")
exec(f"coolant_cell_outer_{code}.temperature = coolant_temp")
exec(f"bundle_universe_{code} = openmc.Universe(cells=(coolant_cell_inner_{code}, coolant_cell_outer_{code}))")
exec(f"pin_universe_{code} = openmc.Universe(cells=(fuel_cell_{code}, clad_side_cell_{code}, clad_back_cell_{code}, clad_front_cell_{code}))")
num_pins = [1, 6, 12, 18]
angles = [0, 0, 15, 0]
ring_radii = np.array([0.0, 1.4885, 2.8755, 4.3305])
for i, (r, n, a) in enumerate(zip(ring_radii, num_pins, angles)):
for j in range(n):
# Determine location of center of pin
theta = (a + j/n*360.) * pi/180.
x = r*cos(theta)
y = r*sin(theta)
pin_boundary = openmc.ZCylinder(x0=x, y0=y, r=r_clad)
eval(f"coolant_cell_inner_{code}").region &= +pin_boundary
eval(f"coolant_cell_outer_{code}").region &= +pin_boundary
# Create each fuel pin -- note that we explicitly assign an ID so
# that we can identify the pin later when looking at tallies
pin = openmc.Cell(fill=eval(f"pin_universe_{code}"), region=-pin_boundary)
pin.translation = (x, y, 0)
eval(f"bundle_universe_{code}").add_cell(pin)
exec(f"bundle_{code} = openmc.Cell(fill=bundle_universe_{code}, region=-pt_inner)")
for i, code in enumerate(fuel_array):
exec(f"pressure_tube_{code} = openmc.Cell(fill=Pressure_tube, region=+pt_inner & -pt_outer & +back_surf_clad & -front_surf_clad)")
exec(f"pressure_tube_{code}.temperature = pressure_tube_temp")
exec(f"CO2_layer_{code} = openmc.Cell(fill=CO2, region =+pt_outer & -calandria_inner & +back_surf_clad & -front_surf_clad)")
exec(f"CO2_layer_{code}.temperature = anulus_temp")
exec(f"calandria_{code} = openmc.Cell(fill=Calandria_tube, region=+calandria_inner & -calandria_outer & +back_surf_clad & -front_surf_clad)")
exec(f"calandria_{code}.temperature = calandria_tube_temp")
for i, code in enumerate(fuel_array):
exec(f"pitch_cell_{code} = openmc.Cell(fill=Moderator, region= +calandria_outer & +left & -right & +bottom & -top & +back & -front)")
for i, code in enumerate(fuel_array):
bundle = eval(f"bundle_{code}")
pressure_tube = eval(f"pressure_tube_{code}")
CO2_layer = eval(f"CO2_layer_{code}")
calandria = eval(f"calandria_{code}")
pitch_cell = eval(f"pitch_cell_{code}")
exec(f"pitch_universe_{code} = openmc.Universe(cells=[bundle, pressure_tube, CO2_layer, calandria, pitch_cell])")
#Pitch
for code in (fuel_array):
Pitch_Universe = eval(f"pitch_universe_{code}")
globals()[code] = Pitch_Universe
lat3d = openmc.RectLattice()
lat3d.lower_left = (-57.15/2, -57.15/2, -24.765)
lat3d.pitch = (pitch, pitch, 49.53)
lat3d.universes = [
[[AAA, BBB],
[CCC, DDD]]]
#defining the reflective boundaries
reflective_distance_x = 28.575
reflective_distance_y = 28.575
void_z = 24.765
x_reflective_plane_left = openmc.XPlane(-reflective_distance_x, boundary_type = 'reflective')
x_reflective_plane_right = openmc.XPlane(+reflective_distance_x, boundary_type = 'reflective')
y_reflective_plane_bottom = openmc.YPlane(-reflective_distance_y, boundary_type = 'reflective')
y_reflective_plane_top = openmc.YPlane(+reflective_distance_y, boundary_type = 'reflective')
z_void_back = openmc.ZPlane(-void_z, boundary_type = 'reflective')
z_void_front = openmc.ZPlane(+void_z, boundary_type = 'reflective')
array_cell = openmc.Cell(fill = lat3d, region = +x_reflective_plane_left & -x_reflective_plane_right & +y_reflective_plane_bottom & -y_reflective_plane_top & +z_void_back & -z_void_front)
root_universe = openmc.Universe(cells=[array_cell])
geometry = openmc.Geometry()
geometry.root_universe = root_universe
geometry.export_to_xml()
bounds = [-reflective_distance_x, -reflective_distance_y, -void_z, reflective_distance_x, reflective_distance_y, void_z]
settings = openmc.Settings()
uniform_dist= openmc.stats.Box(bounds[:3], bounds[3:])
settings.source = openmc.IndependentSource(
space=uniform_dist,
)
settings.batches = batches
settings.inactive = inactive_batches
settings.particles = particles
settings.generations_per_batch = generations
settings.temperature['method'] = 'interpolation'
settings.temperature['tolerance'] = 400
settings.temperature['range'] = [0.0, 2400.0]
settings.output = {'tallies': True}
settings.export_to_xml()
chain = openmc.deplete.Chain.from_xml("/home/smrcenter/openmc/endfb80/chain_casl_pwr.xml")
chain.export_to_xml("chain.xml")
model = openmc.Model(geometry, materials, settings)
return model
crit_ppm, guesses, keffs = openmc.search_for_keff(Boron_build_model, initial_guess=ppm_Boron, model_args=Boron_function_args,
tol=1e-5, print_iterations=True, run_args={'output': False})
print('Critical Boron Concentration: {:.6e}'.format(crit_ppm))
ppm_Boron = crit_ppm
# openmc.reset_auto_ids() # clears ID counters and registries
# if os.path.exists("model.xml"):
# os.remove("model.xml")
materials.export_to_xml()
The depletion runs successfully, the renaming works, but the updated materials are not being used within the criticality search. I believe it’s because I am redefining my materials.xml, so I stopped doing this and didn’t include any material definition within the criticality search but instead it was giving me results like so,
clearly not being able to adjust the boron concentration. Therefore, I decided to remove just the fuel definition as I thought this would be carried forward from the depletion run, and it’s the only part that is changed through depletion so I wouldn’t want to override it, but this gave me the following error,
Furthermore, it doesn’t seem to be updating the boron concentration for the next depletion run although the materials are still taken from the previous depletion. I am probably making a simple mistake, but if someone can help me out I would appreciate it.